Análisis y Predicción Dataset ZILLOW
Rocio Gutiérrez López, Álvaro Caballero Villarreal, Raúl Zornoza Lebrons y Marcos Alberto Pazos Lopez
2020-05-11
1. Resumen ejecutivo
Para llevar a cabo este proyecto hacemos uso del dataset de los precios de venta y la estimación previa de dichos precios en el estado de California, en los condados de Ventura, Orange y Los Angeles, proporcionado por kaggle en la competición llamada Zillow Price en la que lo que se busca predecir es el logerror (logerror=log(Zestimate)−log(SalePrice)).
Por lo que estamos haciendo un análisis residual, la diferencia entre el valor real de la variable dependiente y el valor predicho es lo que se denomina residuo. Todos estos residuos son las porciones sin explicación del modelo actual, por lo que prediciendo las partes que todavía no tienen una explicación correcta mejoraremos el modelo añadiendo las nuevas predicciones encima de las existentes, dando resultado a una predicción mejorada.
En primer lugar, lo que se hace es realizar una limpieza y pre-procesado exhaustiva en la que se examina si el tipo de variable esta correctamente registrado, se crean variables categóricas, se lleva a cabo One Hot Encoding para XGBoost, se tratan los Na´s y se eliminan los outliers. Finalmente, el dataset que se va a utilizar posee más de 100.000 observaciones.
Después de la limpieza intentamos incorporar variables externas para enriquecer nuestro modelo, tales como datos de salud de los condados, indicadores de desempleo, demografía, tasa de crimen, número de viviendas desglosadas en tipo o salario medio de una unidad familiar. Tuvimos varios problemas a la hora de cruzarlos y además nos arrojaban muchos NA´s por lo que finalmente descartamos este proceso.
Para el modelado, en primera instancia la idea era crear un ensamble con varios modelos muy diferentes entre sí entre los cuales estuviera incluido XGBoost, tras hacer varias pruebas se comprobó que los resultados del ensamble eran peores que los del XGBoost de manera aislada. Por lo que finalmente se decidió crear un GBM y un XGBoost que en principio son dos modelos similares, pero nos sirven para llevar acabo una comparativa entre interpretabilidad y grado de ajuste. Estos dos elementos son clave a la hora de seleccionar un modelo de ML y dependiendo del cliente preferirá tener un mejor ajuste, aunque lleve mayor grado de abstracción o poseer un error un poco mayor pero que el modelo se pueda comprender de manera más sencilla.
2. Importación de librerías y carga de datos
library(ggplot2)
library(readxl)
library(readr)
library(yaztheme) #remotes::install_github("joshyazman/yaztheme")
library(leaflet)
library(corrplot)
library(data.table)
library(DT)
require(stringr)
library(tidyverse)
library(caret)
library(tidymodels)
library(skimr)
library(DataExplorer)
library(ggpubr)
library(univariateML)
library(GGally)
library(knitr)
library(dplyr)
library(h2o)
library(DALEX)
library(DALEXtra)
library(univariateML)
library(iml)
library(caret)
library(lubridate)
library(doParallel)
library(xgboost)
library(mltools)
library(data.table)
library(MLmetrics)
library(Matrix)
mipaleta <- c("#ECACF5", "#E47BF3", "#7632FC", "#6377D0", "#A8B8FF")data_dict <- read_excel("./zillow_data_dictionary.xlsx")
all_2016 <- read_csv('./properties_2016.csv')
train_2016 <- read_csv('./train_2016_v2.csv')
all_2017 <- read_csv('./properties_2017.csv')
train_2017 <- read_csv('./train_2017.csv')
full_2016 <- left_join(train_2016, all_2016)
full_2017 <- left_join(train_2017, all_2017)
full_data <- full_join(full_2016, full_2017)2.1 Tabla datasets
Properties 2016
Properties 2017
Dataset que se va a utilizar como base
3. EDA y selección de predictores
!!! - Pulsa en estructura para cargarlo :)
10 primeras filas
Estructura
## Classes 'spec_tbl_df', 'tbl_df', 'tbl' and 'data.frame': 167888 obs. of 60 variables:
## $ parcelid : num 11016594 14366692 12098116 12643413 14432541 ...
## $ logerror : num 0.0276 -0.1684 -0.004 0.0218 -0.005 ...
## $ transactiondate : Date, format: "2016-01-01" "2016-01-01" ...
## $ airconditioningtypeid : num 1 NA 1 1 NA 1 NA NA NA NA ...
## $ architecturalstyletypeid : logi NA NA NA NA NA NA ...
## $ basementsqft : num NA NA NA NA NA NA NA NA NA NA ...
## $ bathroomcnt : num 2 3.5 3 2 2.5 4 1 2.5 1 2 ...
## $ bedroomcnt : num 3 4 2 2 4 4 2 3 2 2 ...
## $ buildingclasstypeid : num NA NA NA NA NA NA NA NA NA NA ...
## $ buildingqualitytypeid : num 4 NA 4 4 NA 1 7 NA NA NA ...
## $ calculatedbathnbr : num 2 3.5 3 2 2.5 4 1 2.5 1 2 ...
## $ decktypeid : num NA NA NA NA NA NA NA NA NA NA ...
## $ finishedfloor1squarefeet : num NA NA NA NA NA NA NA 853 NA NA ...
## $ calculatedfinishedsquarefeet: num 1684 2263 2217 839 2283 ...
## $ finishedsquarefeet12 : num 1684 2263 2217 839 2283 ...
## $ finishedsquarefeet13 : num NA NA NA NA NA NA NA NA NA NA ...
## $ finishedsquarefeet15 : num NA NA NA NA NA NA NA NA NA NA ...
## $ finishedsquarefeet50 : num NA NA NA NA NA NA NA 853 NA NA ...
## $ finishedsquarefeet6 : num NA NA NA NA NA NA NA NA NA NA ...
## $ fips : chr "06037" "06059" "06037" "06037" ...
## $ fireplacecnt : num NA NA NA NA NA NA NA 1 NA NA ...
## $ fullbathcnt : num 2 3 3 2 2 4 1 2 1 2 ...
## $ garagecarcnt : num NA 2 NA NA 2 NA NA 2 1 1 ...
## $ garagetotalsqft : num NA 468 NA NA 598 NA NA 0 0 0 ...
## $ hashottuborspa : logi NA NA NA NA NA NA ...
## $ heatingorsystemtypeid : num 2 NA 2 2 NA 2 7 NA NA NA ...
## $ latitude : num 34280990 33668120 34136312 33755800 33485643 ...
## $ longitude : num -1.18e+08 -1.18e+08 -1.18e+08 -1.18e+08 -1.18e+08 ...
## $ lotsizesquarefeet : num 7528 3643 11423 70859 6000 ...
## $ poolcnt : num NA NA NA NA 1 NA NA NA NA NA ...
## $ poolsizesum : num NA NA NA NA NA NA NA NA NA NA ...
## $ pooltypeid10 : logi NA NA NA NA NA NA ...
## $ pooltypeid2 : logi NA NA NA NA NA NA ...
## $ pooltypeid7 : num NA NA NA NA 1 NA NA NA NA NA ...
## $ propertycountylandusecode : chr "0100" "1" "0100" "010C" ...
## $ propertylandusetypeid : num 261 261 261 266 261 261 261 266 266 266 ...
## $ propertyzoningdesc : chr "LARS" NA "PSR6" "LAR3" ...
## $ rawcensustractandblock : chr "060371066.461001" "060590524.222024" "060374638.003004" "060372963.002002" ...
## $ regionidcity : num 12447 32380 47019 12447 17686 ...
## $ regionidcounty : num 3101 1286 3101 3101 1286 ...
## $ regionidneighborhood : num 31817 NA 275411 54300 NA ...
## $ regionidzip : num 96370 96962 96293 96222 96961 ...
## $ roomcnt : num 0 0 0 0 8 0 0 6 0 5 ...
## $ storytypeid : num NA NA NA NA NA NA NA NA NA NA ...
## $ threequarterbathnbr : num NA 1 NA NA 1 NA NA 1 NA NA ...
## $ typeconstructiontypeid : logi NA NA NA NA NA NA ...
## $ unitcnt : num 1 NA 1 1 NA 1 1 NA NA NA ...
## $ yardbuildingsqft17 : num NA NA NA NA NA NA NA NA NA NA ...
## $ yardbuildingsqft26 : num NA NA NA NA NA NA NA NA NA NA ...
## $ yearbuilt : num 1959 2014 1940 1987 1981 ...
## $ numberofstories : num NA NA NA NA 2 NA NA 2 NA 1 ...
## $ fireplaceflag : logi NA NA NA NA NA NA ...
## $ structuretaxvaluedollarcnt : num 122754 346458 61994 171518 169574 ...
## $ taxvaluedollarcnt : num 360170 585529 119906 244880 434551 ...
## $ assessmentyear : num 2015 2015 2015 2015 2015 ...
## $ landtaxvaluedollarcnt : num 237416 239071 57912 73362 264977 ...
## $ taxamount : num 6736 10153 11484 3049 5489 ...
## $ taxdelinquencyflag : chr NA NA NA NA ...
## $ taxdelinquencyyear : num NA NA NA NA NA NA NA NA NA NA ...
## $ censustractandblock : num 6.04e+13 NA 6.04e+13 6.04e+13 6.06e+13 ...3.1 EDA y correción de variables categóricas
Antes de entrenar un modelo predictivo, es muy importante realizar una exploración descriptiva de los datos que se poseen. Este proceso permite entender mejor qué información contiene cada variable, así como detectar posibles errores y llevar a cabo una limpieza exhaustiva. El objetivo es detectar si hay columnas que están almacenadas con el tipo incorrecto o si hay variables que no tienen sentido, como por ejemplo que en metros cuadrados el piso sea 0, es ilógico.
Cuando se entrena un modelo, es importante incluir como predictores únicamente aquellas variables que están realmente relacionadas con la variable respuesta, ya que son estas las que contienen información útil para el modelo. Incluir un exceso de variables suele conllevar una reducción de la capacidad predictiva del modelo cuando se expone a nuevos datos (overfitting).
Algunos algoritmos de Machine Learning, como por ejemplo Random Forest, Lasso o Boosting, contienen sus propias estrategias para seleccionar predictores. A pesar de ello, se decide llevar a cabo un pre-procesado y una selección de variables, con el fin de eliminar el mayor ruido posible.
Cambio del nombre de las variables, de manera que así se puedan comprender e intuir mucho más rápido.
names(full_data) <- c('id_parcel', 'log_error', 'transaction_date', 'air_conditioning_typeid','architectural_style_typeid',
'basement_sqft', 'bathroom_cnt', 'bedroom_cnt','building_class_typeid','building_quality_typeid',
'calculated_bathroom','deck_typeid','size_finished_living_area_1', 'calculated_finished_living_area',
'finished_living_area','perimeter_living_area','total_area','size_finished_living_area_2',
'base_unfinished_finished_area','federal_information_processing','number_fireplaces',
'number_full_bathrooms','number_garages','number_square_feet_garages','have_hottub_or_spa',
'type_heating_system','latitude','longitude','lot_size_square_feet','number_pools','pool_size_sum',
'pool_type','pool_with_spa_hottub', 'pool_without_spa_hottub','zoning_county_level', 'type_land_use',
'allowed_land_uses','raw_census__and_blockid_1','id_city','id_county','id_neighborhood','id_zip',
'room_cnt', 'story_typeid','number_3_4_bathrooms','type_construction_typeid',
'number_units_structure_built', 'patio_in_yard', 'Storage_shed/building_in_yard',
'year_built','number_stories','fire_place_flag','assessed_value_built_structure',
'total_tax_assessed_value','tax_assessment_year', 'assessed_value_land_area', 'tax_amount',
'tax_delinquency_flag','tax_delinquency_year', 'raw_census__and_blockid_2')
full_data_original <- full_data
full_data_original$latitude = full_data_original$latitude/1000000
full_data_original$longitude = full_data_original$longitude/1000000Revisión de filas duplicadas
## El número de filas duplicadas es: 0Descripción estadística de todas las variables en bruto
Para llevar a cabo el EDA se utiliza el paquete DataExplorer, que para casos como este es muy útil ya que hace que llevar a cabo el EDA y la limpieza sea mucho más fácil que con un summary.
| Name | full_data |
| Number of rows | 167888 |
| Number of columns | 60 |
| _______________________ | |
| Column type frequency: | |
| character | 5 |
| Date | 1 |
| logical | 6 |
| numeric | 48 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| federal_information_processing | 34 | 1.00 | 5 | 5 | 0 | 3 | 0 |
| zoning_county_level | 35 | 1.00 | 1 | 4 | 0 | 90 | 0 |
| allowed_land_uses | 59099 | 0.65 | 2 | 10 | 0 | 2346 | 0 |
| raw_census__and_blockid_1 | 34 | 1.00 | 12 | 16 | 0 | 57853 | 0 |
| tax_delinquency_flag | 163205 | 0.03 | 1 | 1 | 0 | 1 | 0 |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| transaction_date | 0 | 1 | 2016-01-01 | 2017-09-25 | 2016-10-11 | 616 |
Variable type: logical
| skim_variable | n_missing | complete_rate | mean | count |
|---|---|---|---|---|
| architectural_style_typeid | 167888 | 0.00 | NaN | : |
| have_hottub_or_spa | 163984 | 0.02 | 1 | TRU: 3904 |
| pool_type | 166262 | 0.01 | 1 | TRU: 1626 |
| pool_with_spa_hottub | 165610 | 0.01 | 1 | TRU: 2278 |
| type_construction_typeid | 167888 | 0.00 | NaN | : |
| fire_place_flag | 167494 | 0.00 | 1 | TRU: 394 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| id_parcel | 0 | 1.00 | 1.299536e+07 | 3.016071e+06 | 1.071174e+07 | 1.154899e+07 | 1.254060e+07 | 1.421930e+07 | 1.676893e+08 | ▇▁▁▁▁ |
| log_error | 0 | 1.00 | 1.000000e-02 | 1.700000e-01 | -4.660000e+00 | -3.000000e-02 | 1.000000e-02 | 4.000000e-02 | 5.260000e+00 | ▁▁▇▁▁ |
| air_conditioning_typeid | 114100 | 0.32 | 1.810000e+00 | 2.970000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.300000e+01 | ▇▁▁▁▁ |
| basement_sqft | 167795 | 0.00 | 6.953800e+02 | 5.839500e+02 | 3.800000e+01 | 2.800000e+02 | 5.880000e+02 | 8.190000e+02 | 3.560000e+03 | ▇▂▁▁▁ |
| bathroom_cnt | 34 | 1.00 | 2.290000e+00 | 1.000000e+00 | 0.000000e+00 | 2.000000e+00 | 2.000000e+00 | 3.000000e+00 | 2.000000e+01 | ▇▁▁▁▁ |
| bedroom_cnt | 34 | 1.00 | 3.040000e+00 | 1.150000e+00 | 0.000000e+00 | 2.000000e+00 | 3.000000e+00 | 4.000000e+00 | 1.600000e+01 | ▇▃▁▁▁ |
| building_class_typeid | 167857 | 0.00 | 3.970000e+00 | 1.800000e-01 | 3.000000e+00 | 4.000000e+00 | 4.000000e+00 | 4.000000e+00 | 4.000000e+00 | ▁▁▁▁▇ |
| building_quality_typeid | 60715 | 0.64 | 6.020000e+00 | 1.880000e+00 | 1.000000e+00 | 4.000000e+00 | 7.000000e+00 | 7.000000e+00 | 1.200000e+01 | ▁▆▇▃▁ |
| calculated_bathroom | 1832 | 0.99 | 2.310000e+00 | 9.800000e-01 | 1.000000e+00 | 2.000000e+00 | 2.000000e+00 | 3.000000e+00 | 2.000000e+01 | ▇▁▁▁▁ |
| deck_typeid | 166616 | 0.01 | 6.600000e+01 | 0.000000e+00 | 6.600000e+01 | 6.600000e+01 | 6.600000e+01 | 6.600000e+01 | 6.600000e+01 | ▁▁▇▁▁ |
| size_finished_living_area_1 | 154995 | 0.08 | 1.356490e+03 | 6.610700e+02 | 4.400000e+01 | 9.450000e+02 | 1.252000e+03 | 1.615000e+03 | 7.625000e+03 | ▇▃▁▁▁ |
| calculated_finished_living_area | 896 | 0.99 | 1.778630e+03 | 9.403600e+02 | 2.000000e+00 | 1.183000e+03 | 1.541000e+03 | 2.103000e+03 | 3.564000e+04 | ▇▁▁▁▁ |
| finished_living_area | 8369 | 0.95 | 1.752330e+03 | 9.213600e+02 | 2.000000e+00 | 1.172000e+03 | 1.520000e+03 | 2.065000e+03 | 2.192900e+04 | ▇▁▁▁▁ |
| perimeter_living_area | 167813 | 0.00 | 1.395710e+03 | 1.165400e+02 | 1.056000e+03 | 1.344000e+03 | 1.440000e+03 | 1.440000e+03 | 1.584000e+03 | ▁▁▂▇▂ |
| total_area | 161297 | 0.04 | 2.368330e+03 | 1.124300e+03 | 5.600000e+02 | 1.634000e+03 | 2.099000e+03 | 2.845000e+03 | 3.564000e+04 | ▇▁▁▁▁ |
| size_finished_living_area_2 | 154995 | 0.08 | 1.367650e+03 | 6.986700e+02 | 4.400000e+01 | 9.450000e+02 | 1.252000e+03 | 1.619000e+03 | 1.246700e+04 | ▇▁▁▁▁ |
| base_unfinished_finished_area | 167081 | 0.00 | 2.197290e+03 | 1.300550e+03 | 2.570000e+02 | 1.057000e+03 | 1.886000e+03 | 3.240000e+03 | 7.224000e+03 | ▇▆▅▁▁ |
| number_fireplaces | 149992 | 0.11 | 1.190000e+00 | 4.900000e-01 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 5.000000e+00 | ▇▁▁▁▁ |
| number_full_bathrooms | 1832 | 0.99 | 2.250000e+00 | 9.600000e-01 | 1.000000e+00 | 2.000000e+00 | 2.000000e+00 | 3.000000e+00 | 2.000000e+01 | ▇▁▁▁▁ |
| number_garages | 112431 | 0.33 | 1.810000e+00 | 6.000000e-01 | 0.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.400000e+01 | ▇▁▁▁▁ |
| number_square_feet_garages | 112431 | 0.33 | 3.476000e+02 | 2.645700e+02 | 0.000000e+00 | 0.000000e+00 | 4.340000e+02 | 4.880000e+02 | 7.339000e+03 | ▇▁▁▁▁ |
| type_heating_system | 62237 | 0.63 | 3.920000e+00 | 3.640000e+00 | 1.000000e+00 | 2.000000e+00 | 2.000000e+00 | 7.000000e+00 | 2.400000e+01 | ▇▃▁▁▁ |
| latitude | 34 | 1.00 | 3.400678e+07 | 2.651175e+05 | 3.333930e+07 | 3.381292e+07 | 3.402170e+07 | 3.417342e+07 | 3.481877e+07 | ▂▆▇▂▁ |
| longitude | 34 | 1.00 | -1.182011e+08 | 3.600586e+05 | -1.194754e+08 | -1.184130e+08 | -1.181769e+08 | -1.179247e+08 | -1.175546e+08 | ▁▁▆▇▆ |
| lot_size_square_feet | 18442 | 0.89 | 2.951060e+04 | 1.224588e+05 | 1.670000e+02 | 5.702000e+03 | 7.200000e+03 | 1.174675e+04 | 6.971010e+06 | ▇▁▁▁▁ |
| number_pools | 133813 | 0.20 | 1.000000e+00 | 0.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | ▁▁▇▁▁ |
| pool_size_sum | 166050 | 0.01 | 5.189300e+02 | 1.557300e+02 | 2.400000e+01 | 4.200000e+02 | 5.000000e+02 | 6.000000e+02 | 1.750000e+03 | ▁▇▁▁▁ |
| pool_without_spa_hottub | 136112 | 0.19 | 1.000000e+00 | 0.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | ▁▁▇▁▁ |
| type_land_use | 34 | 1.00 | 2.618300e+02 | 5.160000e+00 | 3.100000e+01 | 2.610000e+02 | 2.610000e+02 | 2.660000e+02 | 2.750000e+02 | ▁▁▁▁▇ |
| id_city | 3309 | 0.98 | 3.372513e+04 | 4.692247e+04 | 3.491000e+03 | 1.244700e+04 | 2.521800e+04 | 4.545700e+04 | 3.965560e+05 | ▇▁▁▁▁ |
| id_county | 34 | 1.00 | 2.529620e+03 | 8.037500e+02 | 1.286000e+03 | 1.286000e+03 | 3.101000e+03 | 3.101000e+03 | 3.101000e+03 | ▃▁▁▁▇ |
| id_neighborhood | 100902 | 0.40 | 1.892957e+05 | 1.656877e+05 | 6.952000e+03 | 4.673600e+04 | 1.188720e+05 | 2.747650e+05 | 7.641670e+05 | ▇▆▁▁▁ |
| id_zip | 119 | 1.00 | 9.658639e+04 | 3.723010e+03 | 9.598200e+04 | 9.619300e+04 | 9.638900e+04 | 9.698700e+04 | 3.996750e+05 | ▇▁▁▁▁ |
| room_cnt | 34 | 1.00 | 1.480000e+00 | 2.820000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 1.800000e+01 | ▇▂▁▁▁ |
| story_typeid | 167795 | 0.00 | 7.000000e+00 | 0.000000e+00 | 7.000000e+00 | 7.000000e+00 | 7.000000e+00 | 7.000000e+00 | 7.000000e+00 | ▁▁▇▁▁ |
| number_3_4_bathrooms | 145773 | 0.13 | 1.010000e+00 | 1.100000e-01 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 7.000000e+00 | ▇▁▁▁▁ |
| number_units_structure_built | 58832 | 0.65 | 1.110000e+00 | 9.900000e-01 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 2.370000e+02 | ▇▁▁▁▁ |
| patio_in_yard | 162849 | 0.03 | 3.079200e+02 | 2.274300e+02 | 1.100000e+01 | 1.750000e+02 | 2.520000e+02 | 3.750000e+02 | 3.191000e+03 | ▇▁▁▁▁ |
| Storage_shed/building_in_yard | 167723 | 0.00 | 2.712600e+02 | 2.936100e+02 | 1.200000e+01 | 8.800000e+01 | 1.600000e+02 | 3.200000e+02 | 1.366000e+03 | ▇▂▁▁▁ |
| year_built | 1060 | 0.99 | 1.968570e+03 | 2.378000e+01 | 1.824000e+03 | 1.953000e+03 | 1.970000e+03 | 1.987000e+03 | 2.016000e+03 | ▁▁▂▇▆ |
| number_stories | 129719 | 0.23 | 1.440000e+00 | 5.400000e-01 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 2.000000e+00 | 6.000000e+00 | ▇▁▁▁▁ |
| assessed_value_built_structure | 529 | 1.00 | 1.843460e+05 | 2.192841e+05 | 4.400000e+01 | 8.256800e+04 | 1.340510e+05 | 2.142570e+05 | 1.142179e+07 | ▇▁▁▁▁ |
| total_tax_assessed_value | 36 | 1.00 | 4.726833e+05 | 6.028372e+05 | 2.200000e+01 | 2.026523e+05 | 3.500000e+05 | 5.540020e+05 | 4.906124e+07 | ▇▁▁▁▁ |
| tax_assessment_year | 34 | 1.00 | 2.015460e+03 | 5.000000e-01 | 2.015000e+03 | 2.015000e+03 | 2.015000e+03 | 2.016000e+03 | 2.016000e+03 | ▇▁▁▁▇ |
| assessed_value_land_area | 37 | 1.00 | 2.888805e+05 | 4.456451e+05 | 2.200000e+01 | 8.373000e+04 | 1.977760e+05 | 3.553340e+05 | 4.895220e+07 | ▇▁▁▁▁ |
| tax_amount | 45 | 1.00 | 5.989520e+03 | 7.214730e+03 | 1.992000e+01 | 2.798070e+03 | 4.501040e+03 | 6.914160e+03 | 5.866393e+05 | ▇▁▁▁▁ |
| tax_delinquency_year | 163205 | 0.03 | 1.383000e+01 | 2.420000e+00 | 3.000000e+00 | 1.300000e+01 | 1.400000e+01 | 1.500000e+01 | 9.900000e+01 | ▇▁▁▁▁ |
| raw_census__and_blockid_2 | 886 | 0.99 | 6.049390e+13 | 1.054137e+12 | 6.037101e+13 | 6.037312e+13 | 6.037604e+13 | 6.059042e+13 | 4.830301e+14 | ▇▁▁▁▁ |
Antes de llevar a cabo una gran limpieza de NA´s revisamos las variables categóricas o las que podrían serlo y no estan categorizadas de esta forma, ya que modificándolas un poco podemos hacer que esos NA´s se conviertan en una nivel de la variable factor y así nutrir nuestro modelo. Este proceso tiene especialmente sentido en las variables lógicas.
Análisis de las variables categóricas y lógicas.
plot_bar(
full_data,
ncol = 3, nrow = 4,
ggtheme = theme_minimal(),
theme_config = list(
strip.text = element_text(colour = "#D88EF1", size = 10, face = 2),
legend.position = "none"
)
)
Tax_delinquency_flag
Esta variable significa que los impuestos a la propiedad para este paquete están vencidos a partir de 2015. Esta variable es un claro ejemplo de cómo a veces lo datos vienen mal registrados y si no se examinan con exactitud se puede cometer graves errores. Automáticamente el sistema registra como NA los que no tienen vencido ese impuesto, lo que hay que hacer es convertir los NA en un nivel y dejar una variable que sea dicotómica YES o NO. Además, se aprecia como muy pocas viviendas en la muestra poseen este retraso.
full_data %>%
ggplot(aes(x=as.factor(tax_delinquency_flag), y=(log_error), color = tax_delinquency_flag)) +
geom_jitter(alpha=0.3, color='#99B4E9') +
geom_boxplot(outlier.color="#8E60E3", color='#6800FF') +
labs(x='Impuestos vencidos', title ='Distribución del log error en función de impuestos vencidos')+
theme_minimal() 
Fire_place_flag
full_data <- full_data %>% mutate(fire_place_flag = as.factor(fire_place_flag))
levels <- levels(full_data$fire_place_flag)
levels[length(levels) + 1] <- "FALSE"
full_data$fire_place_flag <- factor(full_data$fire_place_flag, levels = levels)
full_data$fire_place_flag[is.na(full_data$fire_place_flag)] <- "FALSE"Esta variable indica si posee chimenea o no. Por ejemplo esta variable si que tiene sentido binarizarla, es decir, crear nuevas variables dummy con cada uno de los niveles de las variables cualitativas. Pero no se va a llevar a cabo este proceso de manera manual ya que hay modelos que las crean de manera automática en su proceso, como se apreciará más adelante con Xgboost.
Have_hottub_or_spa, pool_type y pool_with_spa/hottub
Las 3 variables coinciden en que tienen muy pocos valores con TRUE y resto con NA, se lleva a cabo el mismo proceso que se llevó con tax_delinquency_flag y los NA´s se convierten en FALSE.
full_data <- full_data %>% mutate(have_hottub_or_spa = as.factor(have_hottub_or_spa))
levels1 <- levels(full_data$have_hottub_or_spa)
levels1[length(levels) + 1] <- "FALSE"
full_data$have_hottub_or_spa <- factor(full_data$have_hottub_or_spa, levels = levels1)
full_data$have_hottub_or_spa[is.na(full_data$have_hottub_or_spa)] <- "FALSE"
full_data <- full_data %>% mutate(pool_type = as.factor(pool_type))
levels2 <- levels(full_data$pool_type)
levels2[length(levels) + 1] <- "FALSE"
full_data$pool_type <- factor(full_data$pool_type, levels = levels2)
full_data$pool_type[is.na(full_data$pool_type)] <- "FALSE"
full_data <- full_data %>% mutate(pool_with_spa_hottub = as.factor(pool_with_spa_hottub))
levels3 <- levels(full_data$pool_with_spa_hottub)
levels3[length(levels) + 1] <- "FALSE"
full_data$pool_with_spa_hottub <- factor(full_data$pool_with_spa_hottub, levels = levels3)
full_data$pool_with_spa_hottub[is.na(full_data$pool_with_spa_hottub)] <- "FALSE"Deck_typeid
Building_class_typeid
full_data <- full_data %>% mutate(building_class_typeid = as.factor(building_class_typeid))
levels11 <- levels(full_data$building_class_typeid)
levels11[length(levels) + 1] <- "ANOTHER"
full_data$building_class_typeid <- factor(full_data$building_class_typeid, levels = levels11)
full_data$building_class_typeid[is.na(full_data$building_class_typeid)] <- "ANOTHER"
skim(full_data$building_class_typeid)| Name | full_data$building_class_… |
| Number of rows | 167888 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| factor | 1 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| data | 0 | 1 | FALSE | 3 | ANO: 167857, 4: 30, 3: 1 |
| Se puede aprecia | r como despu | és convertir la | variable a | categórica | , el nivel en el que más observacioens hay, es del tipo edificios con marcos de madera o madera y acero. Se mantiene esta variable aunque realmente no aporte gran información y la inmensa mayoría están en el nivel another que se ha creado. |
| A pesar de ello | se decide ma | ntenerla, ya que | el modelo | se encarga | rá de desecharla si confirma las sospechas de que realmente no aporta nada. |
Building_quality_typeid
full_data <- full_data %>% mutate(building_quality_typeid = as.factor(building_quality_typeid))
levels12 <- levels(full_data$building_quality_typeid)
levels12[length(levels) + 1] <- "ANOTHER"
full_data$building_quality_typeid <- factor(full_data$building_quality_typeid, levels = levels12)
full_data$building_quality_typeid[is.na(full_data$building_quality_typeid)] <- "ANOTHER"full_data %>%
ggplot(aes(x=as.factor(building_quality_typeid), y=(log_error), color = building_quality_typeid)) +
geom_jitter(alpha=0.3, color='#99B4E9') +
geom_boxplot(outlier.color="#8E60E3", color='#6800FF') +
labs(x='Calidad de construcción', title ='Distribución del log error en función de la calidad de construcción')+
theme_minimal() 
Se observa como es basante intuitivo que los NA´s encajan con el número que falta.
Air_conditioning_typeid
full_data <- full_data %>% mutate(air_conditioning_typeid = as.factor(air_conditioning_typeid))
levels13 <- levels(full_data$air_conditioning_typeid)
levels13[length(levels) + 1] <- "ANOTHER"
full_data$air_conditioning_typeid <- factor(full_data$air_conditioning_typeid, levels = levels13)
full_data$air_conditioning_typeid[is.na(full_data$air_conditioning_typeid)] <- "ANOTHER"full_data %>%
ggplot(aes(x=as.factor(air_conditioning_typeid), y=(log_error), color = air_conditioning_typeid)) +
geom_jitter(alpha=0.3, color='#99B4E9') +
geom_boxplot(outlier.color="#8E60E3", color='#6800FF') +
labs(x='Tipo de aire acondicionado', title ='Distribución del log error en función del tipo de aire acondicionado')+
theme_minimal() 
Number_fireplaces
full_data <- full_data %>% mutate(number_fireplaces = as.factor(number_fireplaces))
levels14 <- levels(full_data$number_fireplaces)
levels14[length(levels) + 1] <- "ANOTHER"
full_data$number_fireplaces <- factor(full_data$number_fireplaces, levels = levels14)
full_data$number_fireplaces[is.na(full_data$number_fireplaces)] <- "ANOTHER"Number_garages
full_data <- full_data %>% mutate(number_garages = as.factor(number_garages))
levels16 <- levels(full_data$number_garages)
levels16[length(levels) + 1] <- "FALSE"
full_data$number_garages <- factor(full_data$number_garages, levels = levels16)
full_data$number_garages[is.na(full_data$number_garages)] <- "FALSE"El número de chimeneas o garajes, son variables que como principal idea obviamente son numéricas, pero la intención es no tratar esos NA´s como un 0 ya que por ejemplo hay 525 viviendas que registraron tener 0 chimeneas, por lo que ese 0 ya esta registrado, si no tratarlas como una variable categórica, en la que esta dentro de una categoría.
Type_heating_system
full_data <- full_data %>% mutate(type_heating_system = as.factor(type_heating_system))
levels18 <- levels(full_data$type_heating_system)
levels18[length(levels) + 1] <- "FALSE"
full_data$type_heating_system <- factor(full_data$type_heating_system, levels = levels18)
full_data$type_heating_system[is.na(full_data$type_heating_system)] <- "FALSE"En este caso, el tipo 2 es calefacción central y la mayoría poseen dicho tipo.
Federal_information_processing
##
## 06037 06059 06111
## 109270 45136 13448Se observa la distribución de los valores de la variable federal_information_processing y se ve como la mayoría de los datos están en una de las 3 categorías. Dado que únicamente existen 34 NA´s se decide eliminar las observaciones y no crear una dummy para esos NA´s.
full_data <- full_data %>% mutate(federal_information_processing = as.factor(federal_information_processing))
full_data = full_data [!is.na(full_data$federal_information_processing),]Se intuye que serán identificadores que corresponderán con cada uno de los condados del dataset.
full_data %>%
ggplot(aes(x=as.factor(federal_information_processing), y=(log_error), color = federal_information_processing)) +
geom_jitter(alpha=0.3, color='#99B4E9') +
geom_boxplot(outlier.color="#8E60E3", color='#6800FF') +
labs(x='Federal information', title ='Distribución del log error en función del federal information')+
theme_minimal() 
Zoning_county_level
## [1] "0100" "1" "010C" "122" "1129" "34" "1128" "010E" "0104" "0101"
## [11] "0200" "0700" "1111" "01DC" "010D" "1110" "0400" "012C" "010V" "1116"
## [21] "01HC" "010G" "0300" "010F" "1117" "0103" "38" "1210" "0111" "010M"
## [31] "96" "135" "0108" "1014" "1112" "0201" "0109" "1310" "010H" "1410"
## [41] "1222" "1321" "1720" "1011" "1432" "0401" "0102" "012D" "73" "105"
## [51] "0110" "100V" "0130" "8800" "0303" "0210" "1012" "1333" "0114" "01DD"
## [61] "020G" "040A" "012E" "020M" "040V" "070D" "1200" "030G" "1722" "6050"
## [71] "1421" "010" NA "200" "0" "1420" "0131" "0301" "01HE" "0204"
## [81] "0113" "0133" "040G" "1120" "0105" "0141" "0115" "010L" "040B" "0203"
## [91] "020E"Esta variable indica los diferentes códigos de uso de la tierra en los condados, que como se puede apreciar existe gran variedad, en concreto hay 90 tipos y 35 Na´s que al ser tan pocos se decide eliminarlos. Se observa como la mayoría de las observaciones pertenecen a los tipos: 0100, 010C, 122. Posee 90 niveles y hay modelos que no van a poder lidiar con ello por lo que se elimina esta variable.
| Name | full_data$zoning_county_l… |
| Number of rows | 167853 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| factor | 1 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| data | 0 | 1 | FALSE | 90 | 010: 57628, 122: 28450, 010: 19102, 010: 13848 |
Allowed_land_uses
Esta variable en principio es similar a la anterior, ya que consiste en la descripción de los usos del suelo permitidos (zonificación) para esa propiedad. Posee 2346 niveles y hay modelos que no van a poder lidiar con ello por lo que se elimina esta variable.
Story_typeid
Patio_in_yard
Esta variable significa tipo de patio que posee, por lo que con los NA´s se crea un tipo adicional llamado another.
full_data <- full_data %>% mutate(patio_in_yard = as.factor(patio_in_yard))
levels9 <- levels(full_data$patio_in_yard)
levels9[length(levels) + 1] <- "ANOTHER"
full_data$patio_in_yard <- factor(full_data$patio_in_yard, levels = levels9)
full_data$patio_in_yard[is.na(full_data$patio_in_yard)] <- "ANOTHER"Number_stories
Número de historias o niveles que tiene la casa, se crea una categoría adicional con los Na´s
full_data <- full_data %>% mutate(number_stories = as.factor(number_stories))
levels10 <- levels(full_data$number_stories)
levels10[length(levels) + 1] <- "ANOTHER"
full_data$number_stories <- factor(full_data$number_stories, levels = levels10)
full_data$number_stories[is.na(full_data$number_stories)] <- "ANOTHER"Después de los NA´s del que más existen es del tipo 1 que es ático y sótano y del tipo 2 que es ático
Tax_assessment_year
Posee 2 niveles o 2015 o 2016.
Raw_census__and_blockid_1 Esta variable que aparece 2 veces y en principio son identicas, para ello cogemos un par de valores aleatorios y observamos si coinciden :raw_census__and_blockid_1 y raw_census__and_blockid_2
Efectivamente coinciden y solamente se diferencian en el formato. Llevando a cabo un análisis profundo, link[https://www.ffiec.gov/census/Default.aspx]. Finalmente se observa como la variable rawcensusandblock esta compuesta por: - FIPS Code - Number - Block Number Por lo que se procede a crear 2 columnas llamadas : _number y _block.
Se borra la columna raw_census__and_blockid_2 y se mantiene la 1 , ya que como posteriormente se analizara la importancia de variables de momento no molesta.
Se convierte a factor
Architectural_style_typeid y type_construction_typeid
Dado que contienen 167.888 NA´s, es decir un 100% se eliminan directamente.
Number_pools y pool_without_spa.hottub
Dado que son un conteo y las pocas casas que tienen solo tienen 1 y el resto son muchísimos NA´s se eliminan, porque esta información ya se ha guardado antes con variables lógicas
Id_parcel
Para predecir no tiene sentido tener en cuenta la variable identificadora de la parcela ya que es un indicador meramente.
Type_land_use
Se detecta que es categórica y posee 14 niveles
full_data <- full_data %>% mutate(type_land_use = as.factor(type_land_use))
unique(full_data$type_land_use)## [1] 261 266 265 246 263 269 248 247 267 275 260 31 264 47
## Levels: 31 47 246 247 248 260 261 263 264 265 266 267 269 275Por ejemplo del tipo que más hay es 261 que es residencia unifamiliar y 266 que es condominio. Se puede llevar a cabo el proceso de binarización en muchas de esta variables, pero se dedice saltarnos ese paso ya que por ejemplo, uno de los modelos que vamos a utilizar, por ser unos de los que más se usan actualmente debido a su buen funcionamiento es XGBoost, este modelo busca automáticamente que variables son categóricas y crea internamente las variables dummy correspondientes.
Creación de columnas año, mes, día y semana.
full_data$year_transaction <- as.factor(substr(as.character(full_data$transaction_date), 1, 4))
full_data$month_transaction <- as.factor(substr(as.character(full_data$transaction_date), 6, 7))
full_data$day_transaction <- as.factor(substr(as.character(full_data$transaction_date), 9, 10))##
## 2016 2017
## 90274 77579##
## 01 02 03 04 05 06 07 08 09 10 11 12
## 13575 12741 17976 17989 20445 22375 19433 20405 14372 4977 1826 1739##
## 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
## 6561 4267 4561 3819 4965 5320 6006 5791 4870 4833 5072 5747 5276 6167 6207 5131
## 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
## 5151 4942 5466 4944 5318 5691 5314 5209 5181 5621 5474 6886 6389 6006 5668## Warning: Unknown or uninitialised column: 'week_transaction'.## < table of extent 0 >3.2 Análisis Missing values
Se decide eliminar las variables que poseen más de un 25% de NA´s ya que consideramos que lo unico que hacen es meterle ruido al modelo y es preferible añadirle variables enriquecidas a través del proceso de ETL.
Junto con el EDA, es básico conocer el número de observaciones disponibles y si todas ellas están completas. Este paso es muy importante ya que algunos algoritmos no aceptan observaciones incompletas o bien se ven muy influenciados por ellas.
nulls <- data.frame(col = as.character(colnames(full_data)),
pct_null = colSums(is.na(full_data))*100/(colSums(is.na(full_data))+colSums(!is.na(full_data))))%>%
filter(col != 'parcelid')
null_count <- ggplot(nulls, aes(x = reorder(col,pct_null), y = pct_null,
fill = ifelse(pct_null > 25, 'Eliminar','Retener')))+
geom_bar(stat = 'identity')+
coord_flip()+
labs(title = 'Distribución de Missing Data',
x = element_blank(), y = 'Porcentaje de Missing Data')+
yaztheme::theme_yaz()+theme(legend.position = 'right')+
geom_hline(yintercept = 25, linetype = 'dashed')+ theme_minimal() +
scale_fill_manual(values = yaz_cols[c(1,6)], name = 'Acción')
null_count
## [1] 167853 443.4 Correlación de variables cuantitativas
Se procede a analizar la correlación existente entre las variables numéricas.
numericas <- select_if(filtrado,is.numeric)
corrplot(cor(numericas, use="na.or.complete"), type="lower", number.digits = 5, col =mipaleta, method = "pie" , rect.col = "pink", tl.col ="#021B8A")
Es importante tener en cuenta que, algunos modelos LM o GLM se ven perjudicados si incorporan predictores altamente correlacionados.
Se observa como las variables que dan datos relacionados con el baño están muy correlacionadas entre sí, lo cual es obvio ya que están aportando la misma información. También, se observa como están totalmente correlacionadas, finished_living_area y calculated_finished_living_area. Además, también lo están tax_amount, total_tax_assessed_value,assessed_value_built_structure y assessed_value_land_area.
Como era previsible el error no posee correlación con el resto de variables.
3.3 Corrección de variables continuas
Se lleva acabo un análisis exploratorio a través de los siguientes gráficos para ver su distribución.
plot_density(
data = filtrado ,
ncol = 3, nrow = 3,
ggtheme = theme_minimal(),
theme_config = list(
strip.text = element_text(colour = "#D88EF1", size = 12, face = 2)
)
)
Una vez llevada a cabo una limpieza bastante grande de las cualitativas se procede a llevar a cabo una limpieza más exhaustiva de outliers de las cuantitativas.
Se observa que latitud y longitud no aparecen en el formato correcto, por lo que se modifican.
3.5 Caracteristicas propias de la vivienda
Se comienza con un análisis de los baños de las viviendas
bathrooms_related <- select (filtrado,c(calculated_bathroom , number_full_bathrooms,bathroom_cnt ))
skim(bathrooms_related)| Name | bathrooms_related |
| Number of rows | 167853 |
| Number of columns | 3 |
| _______________________ | |
| Column type frequency: | |
| numeric | 3 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| calculated_bathroom | 1797 | 0.99 | 2.31 | 0.98 | 1 | 2 | 2 | 3 | 20 | ▇▁▁▁▁ |
| number_full_bathrooms | 1797 | 0.99 | 2.25 | 0.96 | 1 | 2 | 2 | 3 | 20 | ▇▁▁▁▁ |
| bathroom_cnt | 0 | 1.00 | 2.29 | 1.00 | 0 | 2 | 2 | 3 | 20 | ▇▁▁▁▁ |
Recordemos que la variable number_full_bathrooms, considera baños completos los que contienen lavabo, ducha + bañera e inodoro. Además dado que las variables calculated_bathroom y bathroom_cnt contienen la misma información que es el nº de baños que posee la vivienda incluyendo los baños fraccionales, se elimina la variable que mayor nº de NA´s posee, es decir calculated_bathroom.
Visualizamos los gráficos de caja de las 2 variables para detectar outliers y eliminarlos.
Bathroom_cnt
filtrado %>%
ggplot(aes(x=as.factor(bathroom_cnt), y=(log_error), color = bathroom_cnt)) +
geom_jitter(alpha=0.3, color='#99B4E9') +
geom_boxplot(outlier.color="#8E60E3", color='#6800FF') +
labs(x='Número de baños', title='Distribución del log error en función del nº de baños')+
theme_minimal() 
En vista de estos resultados se decide solo tener en cuenta las viviendas que poseen 8 baños o menos ya que para el resto apenas existen datos y se considera que se salen de lo establecido como “normal”.
Number_full_bathrooms
filtrado %>%
ggplot(aes(x=as.factor(number_full_bathrooms), y=(log_error))) +
geom_jitter(alpha=0.3, color='#99B4E9') +
geom_boxplot(outlier.color="#8E60E3", color='#6800FF') +
labs(x='Número de baños completos', title ='Distribución del log error en función del nº de baños')+
theme_minimal() 
Se filtran los que poseen menos de 8 baños.Además, aunque las 2 esten totalmente correlacionadas se decide dejarlas para analizarlas después en las feature importance
Número de habitaciones
filtrado %>%
ggplot(aes(x=as.factor(bedroom_cnt), y=(log_error))) +
geom_jitter(alpha=0.3, color='#99B4E9') +
geom_boxplot(outlier.color="#8E60E3", color='#6800FF') +
labs(x='Número de habtiaciones', title='Distribución del log error en función del nº de habitaciones')+
theme_minimal() 
Se eliminan todas aquellas viviendas que posean más de 9 habitaciones
Número de salas
filtrado %>%
ggplot(aes(x=as.factor(room_cnt), y=(log_error))) +
geom_jitter(alpha=0.3, color='#99B4E9') +
geom_boxplot(outlier.color="#8E60E3", color='#6800FF') +
labs(x='Número de salas', title='Distribución del log error en función del nº de salas')+
theme_minimal() 
Se decide filtrar a partir de 11 habitaciones y se asume que tiene sentido que existan 0 salas, en los loft.
Pies zonas habitables
También se observa como estan totalmente correlacionadas las variables finished_living_area y calculated_finished_living_area.
living_area <- select(filtrado,c(finished_living_area, calculated_finished_living_area))
skim(living_area)| Name | living_area |
| Number of rows | 165734 |
| Number of columns | 2 |
| _______________________ | |
| Column type frequency: | |
| numeric | 2 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| finished_living_area | 6646 | 0.96 | 1745.76 | 892.97 | 2 | 1171 | 1519 | 2062 | 13598 | ▇▁▁▁▁ |
| calculated_finished_living_area | 17 | 1.00 | 1766.02 | 902.15 | 2 | 1181 | 1538 | 2092 | 22741 | ▇▁▁▁▁ |
| Dado que significan lo mismo y la | variable fin | ished_living_are | a posee mu | chos más | NA´s | se deci | de elim | inarla |
Tras hacer un análisis exploratorio a través de internet se estima que de media en EEUU una vivienda posee 2000 feets, por lo que tiene sentido el siguiente resultado.
| Name | filtrado$calculated_finis… |
| Number of rows | 165734 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| data | 17 | 1 | 1766.02 | 902.15 | 2 | 1181 | 1538 | 2092 | 22741 | ▇▁▁▁▁ |
filtrado %>%
ggplot(aes(x=as.factor(calculated_finished_living_area), y=(log_error))) +
geom_jitter(alpha=0.3, color='#99B4E9') +
geom_boxplot(outlier.color="#8E60E3", color='#6800FF') +
labs(x='Número de salas', title='Distribución del log error en función del area habitable')+
theme_minimal() 
Se estima como minimo 200 feets apra una vivienda y como máximo 10.000 feets
A continuación, se analiza la siguiente variable relacionada con el área por metros cuadrados llamada lot_size_square_feet
| Name | filtrado$lot_size_square_… |
| Number of rows | 165676 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| data | 17970 | 0.89 | 29387.35 | 120901 | 167 | 5705 | 7200 | 11736 | 6971010 | ▇▁▁▁▁ |
Esta variable significa pies totales de la parcela dado que se estima que puede existir entorno a un doble en la parcela que en la vivineda en sí se establecen lo limites con concordacia a ello. Además, posee muchos NA´s por lo que se decide eliminar dichas observaciones.
3.6 Variables relacionadas con los impuestos
tax_variables <- select (filtrado,c(total_tax_assessed_value, assessed_value_built_structure,assessed_value_land_area, tax_amount))
skim(tax_variables)| Name | tax_variables |
| Number of rows | 122706 |
| Number of columns | 4 |
| _______________________ | |
| Column type frequency: | |
| numeric | 4 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| total_tax_assessed_value | 2 | 1 | 462691.42 | 500590.83 | 3254.00 | 194804.00 | 355332.00 | 569000.00 | 25381250.0 | ▇▁▁▁▁ |
| assessed_value_built_structure | 108 | 1 | 173150.96 | 173296.17 | 100.00 | 76724.00 | 128922.50 | 210000.00 | 4653658.0 | ▇▁▁▁▁ |
| assessed_value_land_area | 2 | 1 | 289690.04 | 377570.58 | 370.00 | 79694.50 | 205514.00 | 374058.75 | 22335500.0 | ▇▁▁▁▁ |
| tax_amount | 5 | 1 | 5892.53 | 5931.01 | 19.92 | 2812.11 | 4612.62 | 7030.24 | 288524.6 | ▇▁▁▁▁ |
En función de los cuartiles se filtran para evitar outliers y se eliminan los NA´s.
Total_tax_assessed_value
Esta variable es el impuesto total a la propiedad evaluado para ese año de evaluación. Se establece el limite el 1 millón
filtrado = filtrado [!is.na(filtrado$total_tax_assessed_value),]
filtrado <- filtrado %>% filter(total_tax_assessed_value <= 2000000) filtrado %>%
ggplot(aes(x=total_tax_assessed_value)) +
geom_histogram(bins=400, fill="#CF77D2", color='#77BFD2')+
theme_minimal() 
Assessed_value_built_structure
Esta variable es el valor evaluado de la estructura construida en la parcela.
filtrado = filtrado [!is.na(filtrado$assessed_value_built_structure),]
filtrado <- filtrado %>% filter(assessed_value_built_structure <= 500000) filtrado %>%
ggplot(aes(x=assessed_value_built_structure)) +
geom_histogram(bins=400, fill="#CF77D2", color='#77BFD2')+
theme_minimal() 
Assessed_value_land_area
Es similar a la anterior en este caso es el valor evaluado de la superficie de la parcela.
filtrado = filtrado [!is.na(filtrado$assessed_value_land_area),]
filtrado <- filtrado %>% filter(assessed_value_land_area <= 700000)
skim(filtrado$assessed_value_land_area)| Name | filtrado$assessed_value_l… |
| Number of rows | 112165 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| data | 0 | 1 | 215933.7 | 167135.5 | 370 | 68283 | 183982 | 323600 | 7e+05 | ▇▆▃▂▁ |
filtrado %>%
ggplot(aes(x=assessed_value_land_area)) +
geom_histogram(bins=400, fill="#CF77D2", color='#77BFD2')+
theme_minimal() 
3.7 Análisis de la distribución de la variable target: log error
Distribución del error
full_data %>%
ggplot(aes(x=log_error)) +
geom_histogram(bins=400, fill="#DEA1E0")+
theme_bw()+theme(axis.title = element_text(size=16),axis.text = element_text(size=14))+
ylab("Count")+coord_cartesian(x=c(-0.5,0.5))+
labs(title = "Distribución del log error")+
theme_minimal() 
Se aprecia como está infravalorando y sobrevalorando de igual manera, es decir se puede afirmar que la variable target posee una distribución simétrica.
Distribución del log error a lo largo del tiempo
| Name | filtrado$transaction_date |
| Number of rows | 107547 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| Date | 1 |
| ________________________ | |
| Group variables | None |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| data | 0 | 1 | 2016-01-01 | 2017-09-25 | 2016-10-11 | 602 |
Es decir por cada fecha hay varias observaciones y en total existen 602 fechas. Hay que tener en cuenta que para el rango que estábamos utilizando existen un total de 634 días por lo que al eliminar observaciones se asume que hemos borrado datos para 32 días.
Distribución de la media del error a lo largo del tiempo
transactions <- select (filtrado,c(transaction_date, log_error))
transactions[, "transaction_date"] <- as.character(transactions$transaction_date)
transactions$month <- sapply(strsplit(transactions$transaction_date, "-"), function(x) x[2])
transactions %>%
mutate(year_month = make_date(year=year(transaction_date),month=month(transaction_date)) ) %>%
group_by(year_month) %>% summarize(mean_log_error = mean(log_error)) %>%
ggplot(aes(x=year_month,y=mean_log_error)) +
geom_line(size=1, color="#DEA1E0")+geom_point(size=2, color="darkmagenta")+theme_minimal()
En este gráfico se observa como en el último tramo del dataset la tendencia es alcista, se adelanta que como dataset de set se utilizarán los 3 últimos meses de 2017.
Se analiza si existe relación con el número de transacciones llevadas a cabo en esas fechas:
tmp <- transactions %>% mutate(year_month = make_date(year=year(transaction_date),month=month(transaction_date)))
tmp## # A tibble: 107,547 x 4
## transaction_date log_error month year_month
## <chr> <dbl> <chr> <date>
## 1 2016-01-01 0.0276 01 2016-01-01
## 2 2016-01-02 -0.005 01 2016-01-01
## 3 2016-01-02 0.044 01 2016-01-01
## 4 2016-01-03 0.382 01 2016-01-01
## 5 2016-01-03 0.0344 01 2016-01-01
## 6 2016-01-03 -0.045 01 2016-01-01
## 7 2016-01-03 0.002 01 2016-01-01
## 8 2016-01-03 0.044 01 2016-01-01
## 9 2016-01-03 -0.002 01 2016-01-01
## 10 2016-01-03 0.01 01 2016-01-01
## # ... with 107,537 more rowstmp %>%
group_by(year_month) %>% count() %>%
ggplot(aes(x=year_month,y=n)) +
geom_bar(stat="identity", fill="#DEA1E0") +
theme_minimal()
Ya se demostró anteriormente que el log error no posee correlación con ninguna variable, pero con las que sí que poseía un poco de correlación se procede a analizar si pudieran cuadrar sus distribuciones.
custom_corr_plot <- function(variable1, variable2, df, alpha=0.3){
p <- df %>%
mutate(
# Truco para que se ponga el título estilo facet
title = paste(toupper(variable2), "vs", toupper(variable1))
) %>%
ggplot(aes(x = !!sym(variable1), y = !!sym(variable2))) +
geom_point(alpha = alpha) +
# Tendencia no lineal
geom_smooth(se = FALSE, method = "gam", formula = y ~ splines::bs(x, 3)) +
# Tendencia lineal
geom_smooth(se = FALSE, method = "lm", color = "firebrick") +
facet_grid(. ~ title) +
theme_bw() +
theme(strip.text = element_text(colour = "black", size = 10, face = 2),
axis.title = element_blank())
return(p)
}Se seleccionan las variables que estaban un poco más correladas con el log error aunque como vimos ninguna lo estaba.
variables_continuas <- c("bedroom_cnt", "bathroom_cnt")
plots <- map(
.x = variables_continuas,
.f = custom_corr_plot,
variable2 = "log_error",
df = filtrado
)
ggarrange(plotlist = plots, ncol = 2, nrow = 1) %>%
annotate_figure(
top = text_grob("Correlación", face = "bold", size = 12,
x = 0.4)
)
No tiene mucho sentido porque no están nada correladas. Se seleccionan otras 2 relacionadas con el área habitable para demostrar que tampoco lo están
variables_continuas <- c("calculated_finished_living_area", "lot_size_square_feet")
plots <- map(
.x = variables_continuas,
.f = custom_corr_plot,
variable2 = "log_error",
df = filtrado
)
ggarrange(plotlist = plots, ncol = 2, nrow = 1) %>%
annotate_figure(
top = text_grob("Correlación", face = "bold", size = 12,
x = 0.4)
)
3.8 Construcción de pisos
filtrado %>%
ggplot(aes(x=year_built))+geom_line(stat="density", color="#DEA1E0", size=1.2)+theme_minimal()
Se aprecia como la mayoría de las viviendas de este dataset fueron construidas entorno al año 1950
Se analiza como varia el error en función del año en el que se construyo, en primer lugar con respecto al error absoluto
filtrado %>%
group_by(year_built) %>%
summarize(mean_abs_logerror = mean(abs(log_error)),n()) %>%
ggplot(aes(x=year_built,y=mean_abs_logerror))+
geom_smooth(color="grey40")+
geom_point(color="darkorchid1")+coord_cartesian(ylim=c(0,0.25))+theme_bw()
Se ve como las predicciones de viviendas construidas anteriormente a 1910 el error no estima de manera adecuada, en concreto tiende a infravalorar las viviendas, por lo que tal vez seria buena eliminar esas viviendas
Con respecto al error medio:
filtrado %>%
group_by(year_built) %>%
summarize(mean_logerror = mean(log_error),n()) %>%
ggplot(aes(x=year_built,y=mean_logerror))+
geom_smooth(color="grey40")+
geom_point(color="darkorchid1")+coord_cartesian(ylim=c(0,0.25))+theme_bw()
Con respecto al error medio se observa que especialmente son las que se construyeron antes del 1900, y en concreto existe un outlier que se construyó en 1750, por lo que se decide eliminar esas observaciones.
3.9 Análisis de variables relacionadas con la geolocalización
geolocation <- select (filtrado,c(longitude, latitude, id_county, id_city, id_zip, tract_block, tract_number))
skim(geolocation)| Name | geolocation |
| Number of rows | 107319 |
| Number of columns | 7 |
| _______________________ | |
| Column type frequency: | |
| numeric | 7 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| longitude | 0 | 1.00 | -118.20 | 0.35 | -119.48 | -118.38 | -118.16 | -117.96 | -117.55 | ▁▁▅▇▅ |
| latitude | 0 | 1.00 | 34.03 | 0.26 | 33.34 | 33.84 | 34.02 | 34.18 | 34.72 | ▁▆▇▅▂ |
| id_county | 0 | 1.00 | 2579.12 | 777.91 | 1286.00 | 2061.00 | 3101.00 | 3101.00 | 3101.00 | ▃▁▁▁▇ |
| id_city | 1624 | 0.98 | 34783.95 | 52930.08 | 3491.00 | 12447.00 | 24832.00 | 42150.00 | 396556.00 | ▇▁▁▁▁ |
| id_zip | 8 | 1.00 | 96561.00 | 2650.18 | 95982.00 | 96193.00 | 96389.00 | 96993.00 | 399675.00 | ▇▁▁▁▁ |
| tract_block | 0 | 1.00 | 25094.61 | 23499.97 | 0.00 | 4008.00 | 21006.00 | 41002.00 | 95001.00 | ▇▅▂▁▁ |
| tract_number | 0 | 1.00 | 72770.61 | 20542.97 | 10002.00 | 72342.00 | 75331.00 | 79201.10 | 91106.00 | ▁▁▁▆▇ |
## [1] 3101 1286 2061Se observa que existen 3 identificadores de condando: 3101: Los Angeles 1266: Orange 2061: Ventura Lo utilizaremos para segmentar nuestros modelos.
filtrado %>%
ggplot(aes(x=as.factor(id_county), y=(log_error), color = id_county)) +
geom_jitter(alpha=0.3, color='#99B4E9') +
geom_boxplot(outlier.color="#8E60E3", color='#6800FF') +
labs(x='Counties', title ='Distribución del log error en función del condado')+
theme_minimal() 
Se aprecia como el que más dispersión del error posee es el condado de los Estados Unidos.
Se eliminan las observaciones para las que no existe valor en id_city e id_zip
Análisis longitud y latitud
Latitud
tmptrans <- filtrado %>% mutate(overunder = ifelse(log_error<0,"under","over"))
tmptrans %>% ggplot(aes(x=latitude,y=(log_error),color=overunder))+geom_smooth()+theme_bw()+scale_color_brewer(palette="Set3")
Longitud
Se seleccionan las coordendas en las que es error esta sobreestimando y subestimando para dibujar un cuadrado de la localización
Localización en función del error para ver si hay zonas localizadas que poseean un mayor error
tmp <- full_data_original %>% select(id_parcel, longitude, latitude, log_error)
qpal <- colorQuantile("PuRd", tmp$log_error, n = 7)
leaflet(tmp) %>%
addTiles() %>%
addCircleMarkers(stroke=FALSE, color=~qpal(log_error),fillOpacity = 1) %>%
addLegend("bottomright", pal = qpal, values = ~log_error,title = "Log error",opacity = 1) %>%
addMiniMap()3.10 Análisis de variables con varianza próxima a 0
No es conveniente incluir en los modelos variables predictoras que posean una varianza próxima a cero, es decir, predictores que toman solo unos pocos valores, de los cuales, algunos aparecen con muy poca frecuencia. La función nearZeroVar() del paquete caret identifica como predictores potencialmente problemáticos aquellos que tienen un único valor (cero varianza) o que cumplen las dos siguientes condiciones:
Ratio de frecuencias: ratio entre la frecuencia del valor más común y la frecuencia del segundo valor más común. Este ratio tiende a 1 si las frecuencias están equidistribuidas y a valores grandes cuando la frecuencia del valor mayoritario supera por mucho al resto (el denominador es un número decimal pequeño). Valor por defecto freqCut = 95/5.
Porcentaje de valores únicos: número de valores únicos dividido entre el total de muestras (multiplicado por 100). Este porcentaje se aproxima a cero cuanto mayor es la variedad de valores. Valor por defecto uniqueCut = 10.
Columnas a eliminar
## [1] "building_class_typeid" "deck_typeid" "have_hottub_or_spa"
## [4] "pool_type" "pool_with_spa_hottub" "story_typeid"
## [7] "patio_in_yard" "fire_place_flag" "tax_delinquency_flag"Se seleccionan las variables que no están en la lista
3.11 Visualización del dataset limpio
Se visualiza un resumen final de las variables con las que nos hemos quedado
| Name | filtrado |
| Number of rows | 105691 |
| Number of columns | 33 |
| _______________________ | |
| Column type frequency: | |
| Date | 1 |
| factor | 12 |
| numeric | 20 |
| ________________________ | |
| Group variables | None |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| transaction_date | 0 | 1 | 2016-01-01 | 2017-09-25 | 2016-10-12 | 602 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| air_conditioning_typeid | 0 | 1 | FALSE | 5 | ANO: 78972, 1: 23701, 13: 2961, 11: 37 |
| building_quality_typeid | 0 | 1 | FALSE | 12 | ANO: 35071, 7: 27083, 4: 21355, 6: 10299 |
| federal_information_processing | 0 | 1 | FALSE | 3 | 060: 71433, 060: 24991, 061: 9267 |
| number_fireplaces | 0 | 1 | FALSE | 4 | ANO: 92147, 1: 12326, 2: 1214, 4: 4 |
| number_garages | 0 | 1 | FALSE | 9 | FAL: 99481, 1: 5620, 3: 392, 0: 95 |
| type_heating_system | 0 | 1 | FALSE | 10 | FAL: 39696, 2: 39659, 7: 26209, 13: 53 |
| type_land_use | 0 | 1 | FALSE | 13 | 261: 90309, 266: 6640, 246: 3464, 269: 2675 |
| number_stories | 0 | 1 | FALSE | 4 | ANO: 77686, 1: 17593, 2: 10411, 6: 1 |
| tax_assessment_year | 0 | 1 | FALSE | 2 | 201: 56791, 201: 48900 |
| year_transaction | 0 | 1 | FALSE | 2 | 201: 56791, 201: 48900 |
| month_transaction | 0 | 1 | FALSE | 12 | 06: 14114, 08: 12891, 05: 12767, 07: 12296 |
| day_transaction | 0 | 1 | FALSE | 31 | 28: 4470, 01: 4185, 29: 4087, 14: 3837 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| log_error | 0 | 1 | 0.01 | 0.16 | -4.66 | -0.03 | 0.00 | 0.04 | 3.44 | ▁▁▇▁▁ |
| bathroom_cnt | 0 | 1 | 2.14 | 0.80 | 1.00 | 2.00 | 2.00 | 3.00 | 8.00 | ▇▃▁▁▁ |
| bedroom_cnt | 0 | 1 | 3.19 | 0.98 | 0.00 | 3.00 | 3.00 | 4.00 | 9.00 | ▁▇▃▁▁ |
| calculated_finished_living_area | 0 | 1 | 1667.36 | 639.76 | 240.00 | 1216.00 | 1536.00 | 1995.00 | 7612.00 | ▇▅▁▁▁ |
| number_full_bathrooms | 0 | 1 | 2.09 | 0.79 | 1.00 | 2.00 | 2.00 | 3.00 | 8.00 | ▇▂▁▁▁ |
| latitude | 0 | 1 | 34.02 | 0.26 | 33.34 | 33.84 | 34.02 | 34.18 | 34.72 | ▁▆▇▅▂ |
| longitude | 0 | 1 | -118.20 | 0.35 | -119.48 | -118.38 | -118.16 | -117.96 | -117.57 | ▁▁▃▇▅ |
| lot_size_square_feet | 0 | 1 | 7071.89 | 3216.14 | 435.00 | 5325.00 | 6556.00 | 8000.00 | 20000.00 | ▂▇▂▁▁ |
| raw_census__and_blockid_1 | 0 | 1 | 60490082.16 | 212426.99 | 60371011.10 | 60374002.07 | 60375726.00 | 60590421.09 | 61110091.00 | ▇▃▁▁▁ |
| id_city | 0 | 1 | 34783.82 | 52931.04 | 3491.00 | 12447.00 | 24832.00 | 42150.00 | 396556.00 | ▇▁▁▁▁ |
| id_county | 0 | 1 | 2580.65 | 776.54 | 1286.00 | 2061.00 | 3101.00 | 3101.00 | 3101.00 | ▃▁▁▁▇ |
| id_zip | 0 | 1 | 96561.25 | 2670.11 | 95982.00 | 96193.00 | 96387.00 | 96993.00 | 399675.00 | ▇▁▁▁▁ |
| room_cnt | 0 | 1 | 1.79 | 3.02 | 0.00 | 0.00 | 0.00 | 5.00 | 11.00 | ▇▁▁▁▁ |
| year_built | 0 | 1 | 1961.93 | 22.23 | 1900.00 | 1950.00 | 1960.00 | 1978.00 | 2016.00 | ▁▂▇▅▂ |
| assessed_value_built_structure | 0 | 1 | 133605.63 | 84616.37 | 100.00 | 71136.00 | 115900.00 | 177000.00 | 500000.00 | ▇▇▃▁▁ |
| total_tax_assessed_value | 0 | 1 | 337546.83 | 202307.23 | 3254.00 | 173257.00 | 316193.00 | 472000.00 | 1039978.00 | ▇▇▅▂▁ |
| assessed_value_land_area | 0 | 1 | 203941.20 | 156538.85 | 370.00 | 64183.00 | 175070.00 | 307200.00 | 700000.00 | ▇▆▃▂▁ |
| tax_amount | 0 | 1 | 4335.16 | 2284.81 | 19.92 | 2561.42 | 4152.43 | 5886.96 | 9999.94 | ▅▇▇▅▂ |
| tract_number | 0 | 1 | 72647.35 | 20640.94 | 10002.00 | 72322.00 | 75315.00 | 79201.00 | 91106.00 | ▁▁▁▆▇ |
| tract_block | 0 | 1 | 24976.98 | 23409.15 | 0.00 | 4008.00 | 21006.00 | 41001.00 | 95001.00 | ▇▅▂▁▁ |
| Finalmente se guardan los datos li | mpios en csv | para poder gene | rar los modelo | s |
4. Modelos
Se decide crear 2 modelos muy similares entre sí para llevar a cabo una comparativa: GBM y XGBoost.
Tanto XGBoost como GBM siguen el principio del aumento de gradiente. La diferencia principal entre estos dos modelos de predicción es que XGBoost utiliza una formalización del modelo más regularizada para controlar el sobreajuste, lo que le proporciona un mejor rendimiento. Otra de las grandes diferencias entre estos dos modelos es el rendimiento, ya que con XGBoost se han realizado una gran cantidad de mejoras de rendimiento en diferentes partes de la implementación que hacen que haya una gran diferencia entre velocidad y utilización de memoria, ya que usa matrices dispersas con algoritmos de dispersión, estructuras de datos mejoradas para una utilización más eficiente de la memoria caché del procesador, lo que lo hace más rápido; y un mejor soporte multinúcleo que reduce el tiempo de entrenamiento.
4.1 GBM Gradient Boosting Machine
GBM Gradient Boosting Machine
Boosting Machine (BM) es un tipo de modelo obtenido al combinar (ensembling) múltiples modelos sencillos, también conocidos como weak learners. Está combinación se realiza de manera secuencial de tal forma que, cada nuevo modelo que se incorpora al conjunto, intenta corregir los errores de los anteriores modelos. Como resultado de la combinación de múltiples modelos, Boosting Machine consigue aprender relaciones no lineales entre la variable respuesta y los predictores. Si bien los modelos combinados pueden ser muy variados, H2O, al igual que la mayoría de librerías, utiliza como weak learners modelos basados en árboles.
Gradient Boosting Machine (GBM) es una generalización del modelo de Boosting Machine que permite aplicar el método de descenso de gradiente para optimizar cualquier función de coste durante el ajuste del modelo.
El valor predicho por un modelo GBM es la agregación (normalmente la moda en problemas de clasificación y la media en problemas de regresión) de las predicciones de todos los modelos individuales que forman el ensemble.
Dado que el ajuste de los weak learners que forman un GBM se hace de forma secuencial (el árbol i se ajusta a partir de los resultados del árbol i-1), las posibilidades de paralelizar el algoritmo son limitadas. H2O consigue acelerar el proceso en cierta medida paralelizando el ajuste de los árboles individuales.
Utilizamos H2O porque aunque los comandos se ejecuten desde R, los datos se encuentran en el cluster de H2O, no en memoria. Por lo que hace que se solventen los problemas relacionados con la memoria al tener más de 100.000 variables.
Parada temprana: H2O incorpora criterios de parada para que no se continúe mejorando un modelo si ya se ha alcanzado una solución aceptable. Dependiendo del algoritmo, el criterio de parada es distinto.
También hay que remarcar que el proceso de EDA y selección también se podría haber llevado a cabo con H2O pero para el tamaño del dataset generado no ha sido necesario.
http://localhost:54321/flow/index.html
filtrado <- read_csv('Datos_limpios_05_05_2020.csv')
filtrado = subset(filtrado, select = -c(transaction_date) )
filtrado = subset(filtrado, select = -c(X1))
filtrado$air_conditioning_typeid <- as.factor(filtrado$air_conditioning_typeid)
filtrado$building_quality_typeid <- as.factor(filtrado$building_quality_typeid)
filtrado$number_fireplaces <- as.factor(filtrado$number_fireplaces)
filtrado$number_garages <- as.factor(filtrado$number_garages)
filtrado$type_heating_system <- as.factor(filtrado$type_heating_system)
filtrado$number_stories <- as.factor(filtrado$number_stories)
filtrado$month_transaction <- as.factor(filtrado$month_transaction)
filtrado$day_transaction <- as.factor(filtrado$day_transaction)
filtrado$federal_information_processing <- as.factor(filtrado$federal_information_processing)
filtrado$type_land_use <- as.factor(filtrado$type_land_use)
filtrado$day_transaction <- as.factor(filtrado$day_transaction)
filtrado$tax_assessment_year <- as.factor(filtrado$tax_assessment_year)
filtrado$year_transaction <- as.factor(filtrado$year_transaction)
angeles <- filtrado %>% filter(id_county == 3101)
orange_ventura <- filtrado %>% filter(id_county == 1286 | id_county == 2061)
angeles = subset(angeles, select = -c(id_county))
orange_ventura = subset(orange_ventura, select = -c(id_county))Creación de un cluster local con todos los cores disponibles.
h2o.init(ip = "localhost",
# -1 indica que se empleen todos los cores disponibles.
nthreads = -1,
# Máxima memoria disponible para el cluster.
max_mem_size = "32g")##
## H2O is not running yet, starting it now...
##
## Note: In case of errors look at the following log files:
## C:\Users\ALVARO~1\AppData\Local\Temp\RtmpmMdPXc/h2o_Alvaro_Caballero_started_from_r.out
## C:\Users\ALVARO~1\AppData\Local\Temp\RtmpmMdPXc/h2o_Alvaro_Caballero_started_from_r.err
##
##
## Starting H2O JVM and connecting: Connection successful!
##
## R is connected to the H2O cluster:
## H2O cluster uptime: 2 seconds 588 milliseconds
## H2O cluster timezone: Europe/Paris
## H2O data parsing timezone: UTC
## H2O cluster version: 3.28.1.2
## H2O cluster version age: 1 month and 24 days
## H2O cluster name: H2O_started_from_R_Alvaro_Caballero_cer710
## H2O cluster total nodes: 1
## H2O cluster total memory: 28.44 GB
## H2O cluster total cores: 12
## H2O cluster allowed cores: 12
## H2O cluster healthy: TRUE
## H2O Connection ip: localhost
## H2O Connection port: 54321
## H2O Connection proxy: NA
## H2O Internal Security: FALSE
## H2O API Extensions: Amazon S3, Algos, AutoML, Core V3, TargetEncoder, Core V4
## R Version: R version 3.6.2 (2019-12-12)Condandos Ventura y Orange
Se eliminan los datos del cluster por si ya había sido iniciado
h2o.removeAll()
h2o.no_progress()
datos_h2o <- as.h2o(x = orange_ventura, destination_frame = "datos_h2o")La función h2o.splitFrame() realiza particiones aleatorias
set.seed(123)
particiones <- h2o.splitFrame(data = datos_h2o, ratios = c(0.8), seed = 123)
datos_train_h2o <- h2o.assign(data = particiones[[1]], key = "datos_train_h2o")
datos_test_h2o <- h2o.assign(data = particiones[[2]], key = "datos_test_h2o")
datos_train <- as.data.frame(datos_train_h2o)
datos_test <- as.data.frame(datos_test_h2o)Datos balanceados
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.60500 -0.01820 0.00600 0.01579 0.03340 2.46168## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.354000 -0.018200 0.006812 0.017079 0.034343 2.489000Se define la variable respuesta y los predictores.
Optimización de hiperparámetros para el entrenamiento
hiperparametros <- list(
ntrees = c(500, 1000, 2000),
learn_rate = seq(0.01, 0.1, 0.01),
max_depth = seq(2, 15, 1),
sample_rate = seq(0.5, 1.0, 0.2),
col_sample_rate = seq(0.1, 1.0, 0.3))
# Criterios de parada para la búsqueda
search_criteria <- list(
strategy = "RandomDiscrete",
max_models = 50, # Número máximo de combinaciones
max_runtime_secs = 60*10, # Tiempo máximo de búsqueda
stopping_tolerance = 0.001, # Mejora mínima
stopping_rounds = 5,
seed = 123)
grid_gbm <- h2o.grid(
# Algoritmo y parámetros
algorithm = "gbm",
distribution = "gaussian",
# Variable respuesta y predictores
y = var_respuesta,
x = predictores,
# Datos de entrenamiento
training_frame = datos_train_h2o,
# Preprocesado
ignore_const_cols = TRUE,
# Parada temprana
score_tree_interval = 100,
stopping_rounds = 3,
stopping_metric = "rmse",
stopping_tolerance = 0.01,
# Hiperparámetros optimizados
hyper_params = hiperparametros,
# Estrategia de validación para seleccionar el mejor modelo
seed = 123,
nfolds = 3,
# Tipo de búsqueda
search_criteria = search_criteria,
keep_cross_validation_predictions = TRUE,
grid_id = "grid_gbm")Se muestran los modelos ordenados de mayor a menor rmse
resultados_grid_gbm <- h2o.getGrid(
grid_id = "grid_gbm",
sort_by = "rmse",
decreasing = FALSE)
#as.data.frame(resultados_grid_gbm@summary_table)Se reentrena el modelo con los mejores hiperparámetros
mejores_hiperparam <- h2o.getModel(resultados_grid_gbm@model_ids[[1]])@parameters
modelo_gbm <- h2o.gbm(
# Variable respuesta y predictores
y = var_respuesta,
x = predictores,
# Datos de entrenamiento
training_frame = datos_train_h2o,
# Preprocesado
ignore_const_cols = TRUE,
# Hiperparámetros
learn_rate = mejores_hiperparam$learn_rate,
max_depth = mejores_hiperparam$max_depth,
ntrees = mejores_hiperparam$ntrees,
sample_rate = mejores_hiperparam$sample_rate,
# Estrategia de validación (para comparar modelos)
seed = 123,
nfolds = 10,
keep_cross_validation_predictions = TRUE,
model_id = "modelo_gbm")Importancia predictores
## Variable Importances:
## variable relative_importance scaled_importance
## 1 day_transaction 1225.036377 1.000000
## 2 month_transaction 393.279694 0.321035
## 3 calculated_finished_living_area 238.198410 0.194442
## 4 latitude 236.996582 0.193461
## 5 lot_size_square_feet 178.829773 0.145979
## percentage
## 1 0.331537
## 2 0.106435
## 3 0.064465
## 4 0.064139
## 5 0.048397
##
## ---
## variable relative_importance scaled_importance
## 25 type_land_use 8.795951 0.007180
## 26 air_conditioning_typeid 8.743343 0.007137
## 27 tract_number 0.868093 0.000709
## 28 federal_information_processing 0.000000 0.000000
## 29 raw_census__and_blockid_1 0.000000 0.000000
## 30 year_transaction 0.000000 0.000000
## percentage
## 25 0.002380
## 26 0.002366
## 27 0.000235
## 28 0.000000
## 29 0.000000
## 30 0.000000
Gráficos pdp
par(mfrow = c(3, 2))
pdp_plots <- h2o.partialPlot(
object = modelo_gbm,
data = datos_train_h2o,
cols = predictores,
nbins = 31,
plot = TRUE,
plot_stddev = TRUE)
Diagnostico de residuos de entrenamiento
predicciones_train <- h2o.predict(
modelo_gbm,
newdata = datos_train_h2o)
predicciones_train <- h2o.cbind(
datos_train_h2o["log_error"],
predicciones_train)
predicciones_train <- as.data.frame(predicciones_train)
predicciones_train <- predicciones_train %>%
mutate(residuo = predict - log_error)
p1 <- ggplot(predicciones_train, aes(x = log_error, y = predict)) +
geom_point(alpha = 0.1) +
geom_smooth(method = "gam", color = "red", size = 1, se = FALSE) +
labs(title = "Predicciones vs valor real") +
theme_bw()
p2 <- ggplot(predicciones_train, aes(1:nrow(predicciones_train), y = residuo)) +
geom_point(alpha = 0.1) +
geom_hline(yintercept = 0, color = "red", size = 1) +
labs(title = "Residuos del modelo") +
theme_bw()
p3 <- ggplot(predicciones_train, aes(x = residuo)) +
geom_density() +
geom_rug() +
labs(title = "Residuos del modelo") +
theme_bw()
p4 <- ggplot(predicciones_train, aes(sample = predict)) +
stat_qq() +
stat_qq_line(color = "red", size = 1) +
labs(title = "QQ-plot residuos del modelo") +
theme_bw()
ggpubr::ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2) %>%
ggpubr::annotate_figure(
top = ggpubr::text_grob("Diagnóstico residuos entrenamiento",
color = "black", face = "bold", size = 14))
Predicción
predicciones_test <- h2o.predict(object = modelo_gbm, newdata = datos_test_h2o)
predicciones_test <- as.vector(predicciones_test)
datos_test$prediccion <- predicciones_test
rmse_test_gbm <- MLmetrics::RMSE( y_pred = datos_test$log_error, y_true = datos_test$prediccion)
mse_test_gbm <- MLmetrics::MSE( y_pred = datos_test$log_error, y_true = datos_test$prediccion)
MAE_test_gbm <- MLmetrics::MAE(y_pred = datos_test$log_error,y_true = datos_test$prediccion)
MedianAE_test_gbm <- MLmetrics::MedianAE(y_pred = datos_test$log_error,y_true = datos_test$prediccion)## [1] "Error de test (RMSE) del modelo GBM: 0.160079833705384"## [1] "Error de test (MSE) del modelo GBM: 0.0256255531591432"## [1] "Error de test (MAE) del modelo GBM: 0.0601265360956686"## [1] "Error de test (MedianAE) del modelo GBM: 0.0291977630654353"Condando Los Angeles
Se eliminan los datos del cluster por si ya había sido iniciado
La función h2o.splitFrame() realiza particiones aleatorias
set.seed(123)
particiones <- h2o.splitFrame(data = datos_h2o, ratios = c(0.8), seed = 123)
datos_train_h2o <- h2o.assign(data = particiones[[1]], key = "datos_train_h2o")
datos_test_h2o <- h2o.assign(data = particiones[[2]], key = "datos_test_h2o")
datos_train <- as.data.frame(datos_train_h2o)
datos_test <- as.data.frame(datos_test_h2o)Datos balanceados
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.605000 -0.031498 0.004369 0.012507 0.040200 3.443000## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.655420 -0.031500 0.004372 0.012779 0.039200 3.436000Se define la variable respuesta y los predictores.
Optimización de hiperparámetros para el entrenamiento
hiperparametros <- list(
ntrees = c(500, 1000, 2000),
learn_rate = seq(0.01, 0.1, 0.01),
max_depth = seq(2, 15, 1),
sample_rate = seq(0.5, 1.0, 0.2),
col_sample_rate = seq(0.1, 1.0, 0.3))
# Criterios de parada para la búsqueda
search_criteria <- list(
strategy = "RandomDiscrete",
max_models = 50, # Número máximo de combinaciones
max_runtime_secs = 60*10, # Tiempo máximo de búsqueda
stopping_tolerance = 0.001, # Mejora mínima
stopping_rounds = 5,
seed = 123)
grid_gbm <- h2o.grid(
# Algoritmo y parámetros
algorithm = "gbm",
distribution = "gaussian",
# Variable respuesta y predictores
y = var_respuesta,
x = predictores,
# Datos de entrenamiento
training_frame = datos_train_h2o,
# Preprocesado
ignore_const_cols = TRUE,
# Parada temprana
score_tree_interval = 100,
stopping_rounds = 3,
stopping_metric = "rmse",
stopping_tolerance = 0.01,
# Hiperparámetros optimizados
hyper_params = hiperparametros,
# Estrategia de validación para seleccionar el mejor modelo
seed = 123,
nfolds = 3,
# Tipo de búsqueda
search_criteria = search_criteria,
keep_cross_validation_predictions = TRUE,
grid_id = "grid_gbm")Se muestran los modelos ordenados de mayor a menor rmse
resultados_grid_gbm <- h2o.getGrid(
grid_id = "grid_gbm",
sort_by = "rmse",
decreasing = FALSE)
#as.data.frame(resultados_grid_gbm@summary_table)Se reentrena el modelo con los mejores hiperparámetros
mejores_hiperparam <- h2o.getModel(resultados_grid_gbm@model_ids[[1]])@parameters
modelo_gbm <- h2o.gbm(
# Variable respuesta y predictores
y = var_respuesta,
x = predictores,
# Datos de entrenamiento
training_frame = datos_train_h2o,
# Preprocesado
ignore_const_cols = TRUE,
# Hiperparámetros
learn_rate = mejores_hiperparam$learn_rate,
max_depth = mejores_hiperparam$max_depth,
ntrees = mejores_hiperparam$ntrees,
sample_rate = mejores_hiperparam$sample_rate,
# Estrategia de validación (para comparar modelos)
seed = 123,
nfolds = 10,
keep_cross_validation_predictions = TRUE,
model_id = "modelo_gbm")Importancia predictores
## Variable Importances:
## variable relative_importance scaled_importance
## 1 day_transaction 1921.879395 1.000000
## 2 tax_amount 470.429657 0.244776
## 3 month_transaction 460.579956 0.239651
## 4 calculated_finished_living_area 420.000122 0.218536
## 5 lot_size_square_feet 322.234894 0.167667
## percentage
## 1 0.328925
## 2 0.080513
## 3 0.078827
## 4 0.071882
## 5 0.055150
##
## ---
## variable relative_importance scaled_importance percentage
## 22 type_heating_system 35.815292 0.018636 0.006130
## 23 tax_assessment_year 5.779536 0.003007 0.000989
## 24 air_conditioning_typeid 2.528981 0.001316 0.000433
## 25 number_full_bathrooms 0.000000 0.000000 0.000000
## 26 number_stories 0.000000 0.000000 0.000000
## 27 year_transaction 0.000000 0.000000 0.000000
Gráficos pdp
par(mfrow = c(3, 2))
pdp_plots <- h2o.partialPlot(
object = modelo_gbm,
data = datos_train_h2o,
cols = predictores,
nbins = 31,
plot = TRUE,
plot_stddev = TRUE)
Diagnostico de residuos de entrenamiento
predicciones_train <- h2o.predict(
modelo_gbm,
newdata = datos_train_h2o)
predicciones_train <- h2o.cbind(
datos_train_h2o["log_error"],
predicciones_train)
predicciones_train <- as.data.frame(predicciones_train)
predicciones_train <- predicciones_train %>%
mutate(residuo = predict - log_error)
p1 <- ggplot(predicciones_train, aes(x = log_error, y = predict)) +
geom_point(alpha = 0.1) +
geom_smooth(method = "gam", color = "red", size = 1, se = FALSE) +
labs(title = "Predicciones vs valor real") +
theme_bw()
p2 <- ggplot(predicciones_train, aes(1:nrow(predicciones_train), y = residuo)) +
geom_point(alpha = 0.1) +
geom_hline(yintercept = 0, color = "red", size = 1) +
labs(title = "Residuos del modelo") +
theme_bw()
p3 <- ggplot(predicciones_train, aes(x = residuo)) +
geom_density() +
geom_rug() +
labs(title = "Residuos del modelo") +
theme_bw()
p4 <- ggplot(predicciones_train, aes(sample = predict)) +
stat_qq() +
stat_qq_line(color = "red", size = 1) +
labs(title = "QQ-plot residuos del modelo") +
theme_bw()
ggpubr::ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2) %>%
ggpubr::annotate_figure(
top = ggpubr::text_grob("Diagnóstico residuos entrenamiento",
color = "black", face = "bold", size = 14))
predicciones_test <- h2o.predict(object = modelo_gbm, newdata = datos_test_h2o)
predicciones_test <- as.vector(predicciones_test)
datos_test$prediccion <- predicciones_test
rmse_test_gbm <- MLmetrics::RMSE( y_pred = datos_test$log_error, y_true = datos_test$prediccion)
mse_test_gbm <- MLmetrics::MSE( y_pred = datos_test$log_error, y_true = datos_test$prediccion)
MAE_test_gbm <- MLmetrics::MAE(y_pred = datos_test$log_error,y_true = datos_test$prediccion)
MAPE_test_gbm <- MLmetrics::MAPE(y_pred = datos_test$log_error,y_true = datos_test$prediccion)
MedianAE_test_gbm <- MLmetrics::MedianAE(y_pred = datos_test$log_error,y_true = datos_test$prediccion)## [1] "Error de test (RMSE) del modelo GBM: 0.166860079040484"## [1] "Error de test (MSE) del modelo GBM: 0.0278422859773965"## [1] "Error de test (MAE) del modelo GBM: 0.0741411596768036"## [1] "Error de test (MAPE) del modelo GBM: 22.4478696012399"## [1] "Error de test (MedianAE) del modelo GBM: 0.0372846250869881"4.2 XGBoost
Xgboost significa Extreme Gradient Boosting; este es una implementación específica del método Gradient Boosting Machine, que utiliza aproximaciones más precisas para encontrar el mejor modelo de árbol. Emplea una serie de algoritmos (trucos) que lo hacen realmente exitoso, sobre todo para datos estructurados. Los más importantes son los siguientes:
Calcula gradientes de segundo orden, es decir, segundas derivadas parciales de la pérdida (parecido al método de Newton), lo que nos proporciona información de la dirección de los gradientes y la forma de llegar al mínimo de la función de pérdida. Si bien el aumento de gradiente regular utiliza la función de pérdida de un modelo base (como, por ejemplo, árbol de decisión), como intermediario para minimizar el error del modelo general, XGBoost usa la segunda derivada como una aproximación.
XGBoost tiene una regularización avanzada que mejora la generación del modelo. Como se describe previamente, XGBoost tiene otras ventajas, que el entrenamiento es muy rápido, se puede paralelizar/distribuir a través de clusters.
filtrado$year_built <- as.factor(filtrado$year_built)
filtrado$week_transaction <- as.factor(filtrado$week_transaction)
datos_limpios <- filtradoSacamos las clases de las variables.
feature_classes <- sapply(names(datos_limpios), function(x) {
class(datos_limpios[[x]])})
#feature_classesVariables numéricas.
## [1] "log_error" "bathroom_cnt"
## [3] "bedroom_cnt" "calculated_finished_living_area"
## [5] "number_full_bathrooms" "latitude"
## [7] "longitude" "lot_size_square_feet"
## [9] "raw_census__and_blockid_1" "id_city"
## [11] "id_county" "id_zip"
## [13] "room_cnt" "assessed_value_built_structure"
## [15] "total_tax_assessed_value" "assessed_value_land_area"
## [17] "tax_amount" "tract_number"
## [19] "tract_block"Seleccionamos las variables categóricas.
## [1] "air_conditioning_typeid" "building_quality_typeid"
## [3] "federal_information_processing" "number_fireplaces"
## [5] "number_garages" "type_heating_system"
## [7] "type_land_use" "year_built"
## [9] "number_stories" "tax_assessment_year"
## [11] "year_transaction" "month_transaction"
## [13] "day_transaction" "week_transaction"Creamos una tabla con las variables numericas.
One hot Enconding (Variables dummy para las variables categóricas)
total_categorical <- datos_limpios[, categorical_features]
dummy <- dummyVars(" ~ .", data = total_categorical)
#dataframe con las variables dummy
total_categoric <- data.frame(predict(dummy, newdata = total_categorical))Se crean n total de 277 variables categóricas
Juntamos variables numéricas y las que hemos creado Dummy
Para crear XGboost poseemos un total de 296 variables
Creamos dos datasets diferenciados por su County, uno para Los Angeles donde hay más observaciones, y luego uno formado por Orange y Ventura counties, que además de tener menos obervaciones se diferencian en el median family income, siendo mayor que LA, y pudiendo resultar en un modelo mejor.
datos_LA <- total %>% filter(id_county == 3101)
datos_or_ve <- total %>% filter(id_county == 1286 | id_county == 2061)Creamos ya las particiones train y test, un 80% para train y un 20% para test.
set.seed(123)
spec = c(train = .8, test = .2)
g = sample(cut(
seq(nrow(total)),
nrow(datos_or_ve)*cumsum(c(0,spec)),
labels = names(spec)))
res_orve = split(datos_or_ve, g)
g = sample(cut(
seq(nrow(total)),
nrow(datos_LA)*cumsum(c(0,spec)),
labels = names(spec)))
res_la = split(datos_LA, g)Creamos una variable target y luego el conjunto de variables predictoras.
Creamos una variable Outcome, que tendrá los datos del target.
Creamos accesos directos a los dataset train y test, y les asignamos la clase matriz para introducirla a la matriz xgboost.
total_train_la <- res_la$train
total_test_la <- res_la$test
total_train_orve <- res_orve$train
total_test_orve <- res_orve$testSparse Matrix
train_sparse_la <- sparse.model.matrix(~. , data = total_train_la)[,-1]
test_sparse_la <- sparse.model.matrix(~. , data = total_test_la)[,-1]
train_sparse_orve <- sparse.model.matrix(~. , data = total_train_orve)[,-1]
test_sparse_orve <- sparse.model.matrix(~. , data = total_test_orve)[,-1]Creamos las matrices de XGBoost
dtrain_la <- xgb.DMatrix(train_sparse_la, label = outcome_la)
dtest_la <- xgb.DMatrix(test_sparse_la)
dtrain_orve <- xgb.DMatrix(train_sparse_orve, label = outcome_orve)
dtest_orve <- xgb.DMatrix(test_sparse_orve)
set.seed(1235)Primero entrenamos un Modelo XGBoost Cross Validation que nos permite mejorar de manera notoria los resultados de predicción gracias a la regularización que permite este modelo y que le diferencia de GBM, además de ofrecer el mejor boosting de gradiente. En este primer modelo se realizará la cross-validation, y en el siguiente sin ella, para comprobar si los buenos resultados de la primera son demasiado optimistas.
Primero para Orange y Ventura counties.
param <- list(booster = "gblinear",
objective = "reg:linear",
eval_metric = "mae",
eta=0.5,
gamma = 0.01,
alpha = 0.9,
lambda = 0.8)
set.seed(123)
modelo_orve <- xgb.cv( params = param, data = dtrain_orve, nrounds = 1500, nfold = 5, showsd = T, stratified = T, print_every_n = 500, early_stop_round = 50, maximize = F, prediction = T)## [1] train-mae:0.069867+0.000847 test-mae:0.069902+0.003193
## [501] train-mae:0.056162+0.000870 test-mae:0.056303+0.002911
## [1001] train-mae:0.056162+0.000870 test-mae:0.056303+0.002911
## [1500] train-mae:0.056162+0.000870 test-mae:0.056303+0.002911En este caso, sacamos con el modelo optimizado un error absoluto medio de 0.05 para la validación cruzada en los counties de Orange y Ventura.
## [1] 0.0563024Y ahora para Los Angeles County
param <- list(booster = "gblinear",
objective = "reg:linear",
eval_metric = "mae",
eta=0.5,
gamma = 0.01,
alpha = 0.9,
lambda = 0.8)
set.seed(123)
modelo_la <- xgb.cv( params = param, data = dtrain_la, nrounds = 1500, nfold = 5, showsd = T, stratified = T, print_every_n = 500, early_stop_round = 50, maximize = F, prediction = T)## [1] train-mae:0.086781+0.000242 test-mae:0.086781+0.000948
## [501] train-mae:0.073676+0.000302 test-mae:0.073692+0.000860
## [1001] train-mae:0.073671+0.000302 test-mae:0.073687+0.000859
## [1500] train-mae:0.073669+0.000302 test-mae:0.073686+0.000860En este caso, sacamos con el modelo optimizado un error absoluto medio de 0.07 para la validación cruzada en Los Angeles.
## [1] 0.0736412Vamos a tratar de mejorar el modelo de Validación cruzada de XGBoost a través de una búsqueda de metaparámetros con GridSearch. Creamos un Grid con los metaparámetros para los modelos de OR y VEN, como para el de LA.
searchGridSubCol <- expand.grid(gamma = c(0,0.01,0.1),
eta = c(0.1, 0.3, 0.5),
alpha = c(0.7,0.8,0.9))Grid Search para dataset de Orange y Ventura.
set.seed(123)
system.time(
rmseErrorsHyperparameters <- apply(searchGridSubCol, 1, function(parameterList){
#Extraemos los parámetros para testear
currentEta <- parameterList[["eta"]]
currentGamma <- parameterList[["gamma"]]
currentAlpha <- parameterList[["alpha"]]
param = list(eta = currentEta,
gamma = currentGamma,
alpha = currentAlpha,
eval_metric = "mae",
objective = "reg:linear",
booster = "gblinear")
#hacemos un modelo XGBoost Cross Validation
xgboostModelCV <- xgb.cv(data = dtrain_orve, params = param, nrounds = 5000,
nfold = 5, showsd = TRUE,
metrics = "mae", verbose = T,
print_every_n = 1000, early_stopping_rounds = 100)
xvalidationScores <- as.data.frame(xgboostModelCV$evaluation_log)
mae <- tail(xvalidationScores$test_mae_mean, 1)
trmae <- tail(xvalidationScores$train_mae_mean,1)
stop_it <- tail(xvalidationScores$st)
output <- return(c(mae, trmae, stop_it, currentEta, currentGamma, currentAlpha))}))## [1] train-mae:0.139408+0.000972 test-mae:0.139440+0.002490
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## [1001] train-mae:0.056196+0.000873 test-mae:0.056336+0.002916
## Stopping. Best iteration:
## [923] train-mae:0.056197+0.000873 test-mae:0.056335+0.002917
##
## [1] train-mae:0.139402+0.000533 test-mae:0.139414+0.003349
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## [1001] train-mae:0.056220+0.000497 test-mae:0.056335+0.002353
## Stopping. Best iteration:
## [1320] train-mae:0.056218+0.000497 test-mae:0.056332+0.002354
##
## [1] train-mae:0.139404+0.000747 test-mae:0.139413+0.002844
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## [1001] train-mae:0.056200+0.000741 test-mae:0.056239+0.003315
## Stopping. Best iteration:
## [1334] train-mae:0.056198+0.000741 test-mae:0.056237+0.003314
##
## [1] train-mae:0.073223+0.000602 test-mae:0.073253+0.002060
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [489] train-mae:0.056201+0.000691 test-mae:0.056258+0.002123
##
## [1] train-mae:0.073221+0.001005 test-mae:0.073217+0.003364
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [435] train-mae:0.056191+0.000988 test-mae:0.056246+0.003076
##
## [1] train-mae:0.073219+0.000534 test-mae:0.073245+0.002454
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [288] train-mae:0.056205+0.000581 test-mae:0.056299+0.002611
##
## [1] train-mae:0.069654+0.000998 test-mae:0.069738+0.003265
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [141] train-mae:0.056221+0.001174 test-mae:0.056328+0.003868
##
## [1] train-mae:0.069653+0.000778 test-mae:0.069706+0.002573
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [165] train-mae:0.056208+0.000886 test-mae:0.056321+0.002575
##
## [1] train-mae:0.069659+0.000972 test-mae:0.069698+0.003172
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [187] train-mae:0.056231+0.000895 test-mae:0.056288+0.003011
##
## [1] train-mae:0.142308+0.000973 test-mae:0.142326+0.003557
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## [1001] train-mae:0.056235+0.001000 test-mae:0.056317+0.004058
## Stopping. Best iteration:
## [948] train-mae:0.056235+0.001001 test-mae:0.056317+0.004057
##
## [1] train-mae:0.142309+0.001330 test-mae:0.142346+0.002490
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## [1001] train-mae:0.056226+0.001202 test-mae:0.056358+0.004373
## Stopping. Best iteration:
## [1186] train-mae:0.056224+0.001202 test-mae:0.056357+0.004374
##
## [1] train-mae:0.142305+0.000912 test-mae:0.142303+0.003941
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## [1001] train-mae:0.056212+0.000723 test-mae:0.056276+0.003111
## Stopping. Best iteration:
## [1032] train-mae:0.056212+0.000723 test-mae:0.056276+0.003111
##
## [1] train-mae:0.072978+0.000134 test-mae:0.072996+0.000910
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [367] train-mae:0.056228+0.000253 test-mae:0.056350+0.000979
##
## [1] train-mae:0.072979+0.000623 test-mae:0.073002+0.002604
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [355] train-mae:0.056246+0.000613 test-mae:0.056321+0.002120
##
## [1] train-mae:0.072978+0.000330 test-mae:0.072971+0.001457
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [431] train-mae:0.056205+0.000346 test-mae:0.056244+0.001571
##
## [1] train-mae:0.069800+0.000813 test-mae:0.069812+0.003113
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [185] train-mae:0.056202+0.000736 test-mae:0.056239+0.002361
##
## [1] train-mae:0.069806+0.000681 test-mae:0.069836+0.002428
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [181] train-mae:0.056219+0.000884 test-mae:0.056286+0.002602
##
## [1] train-mae:0.069802+0.001051 test-mae:0.069829+0.003993
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [170] train-mae:0.056221+0.001218 test-mae:0.056310+0.004787
##
## [1] train-mae:0.144666+0.001476 test-mae:0.144741+0.004141
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## [1001] train-mae:0.056247+0.001380 test-mae:0.056364+0.004901
## Stopping. Best iteration:
## [1420] train-mae:0.056245+0.001380 test-mae:0.056361+0.004901
##
## [1] train-mae:0.144659+0.000769 test-mae:0.144678+0.002982
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## [1001] train-mae:0.056237+0.000836 test-mae:0.056346+0.003140
## Stopping. Best iteration:
## [1214] train-mae:0.056235+0.000836 test-mae:0.056344+0.003144
##
## [1] train-mae:0.144659+0.001591 test-mae:0.144643+0.005433
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## [1001] train-mae:0.056230+0.001081 test-mae:0.056285+0.003830
## Stopping. Best iteration:
## [1430] train-mae:0.056227+0.001081 test-mae:0.056283+0.003831
##
## [1] train-mae:0.072921+0.000418 test-mae:0.072951+0.001469
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [358] train-mae:0.056232+0.000337 test-mae:0.056395+0.001205
##
## [1] train-mae:0.072920+0.000652 test-mae:0.072954+0.002427
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [414] train-mae:0.056259+0.000755 test-mae:0.056327+0.002775
##
## [1] train-mae:0.072919+0.000446 test-mae:0.072934+0.001785
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [309] train-mae:0.056223+0.000303 test-mae:0.056325+0.001111
##
## [1] train-mae:0.069864+0.000481 test-mae:0.069901+0.001664
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [384] train-mae:0.056177+0.000436 test-mae:0.056232+0.001764
##
## [1] train-mae:0.069865+0.000547 test-mae:0.069906+0.002137
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [375] train-mae:0.056185+0.000695 test-mae:0.056276+0.002269
##
## [1] train-mae:0.069865+0.000676 test-mae:0.069888+0.002366
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [353] train-mae:0.056152+0.000925 test-mae:0.056287+0.002567## user system elapsed
## 406.34 190.23 240.51Y vemos los resultados ordenados de forma creciente por TestMAE para ver los metaparámetros que dieron un MAE más bajo en el test
set.seed(123)
output <- as.data.frame(t(rmseErrorsHyperparameters))
varnames <- c("TestMAE", "TrainMAE", "Eta", "Gamma", "Alpha")
names(output) <- varnames
output[order(output$TestMAE),]## TestMAE TrainMAE Eta Gamma Alpha
## 25 0.0562332 0.0561776 0.5 0.00 0.9
## 3 0.0562376 0.0561974 0.1 0.10 0.7
## 16 0.0562404 0.0562022 0.5 0.00 0.8
## 15 0.0562448 0.0562054 0.3 0.10 0.8
## 5 0.0562466 0.0561908 0.3 0.01 0.7
## 4 0.0562582 0.0562010 0.3 0.00 0.7
## 12 0.0562760 0.0562110 0.1 0.10 0.8
## 26 0.0562768 0.0561846 0.5 0.01 0.9
## 21 0.0562828 0.0562270 0.1 0.10 0.9
## 17 0.0562862 0.0562194 0.5 0.01 0.8
## 27 0.0562880 0.0561520 0.5 0.10 0.9
## 9 0.0562888 0.0562318 0.5 0.10 0.7
## 6 0.0562998 0.0562020 0.3 0.10 0.7
## 18 0.0563112 0.0562212 0.5 0.10 0.8
## 10 0.0563166 0.0562340 0.1 0.00 0.8
## 14 0.0563224 0.0562450 0.3 0.01 0.8
## 8 0.0563246 0.0562096 0.5 0.01 0.7
## 24 0.0563256 0.0562212 0.3 0.10 0.9
## 23 0.0563278 0.0562590 0.3 0.01 0.9
## 7 0.0563300 0.0562190 0.5 0.00 0.7
## 2 0.0563320 0.0562180 0.1 0.01 0.7
## 1 0.0563356 0.0561958 0.1 0.00 0.7
## 20 0.0563446 0.0562348 0.1 0.01 0.9
## 13 0.0563506 0.0562280 0.3 0.00 0.8
## 11 0.0563568 0.0562240 0.1 0.01 0.8
## 19 0.0563616 0.0562452 0.1 0.00 0.9
## 22 0.0563956 0.0562312 0.3 0.00 0.9Ahora probamos el mismo GridSearch para Los Angeles.
#Creamos un Grid para los modelos de OR y VEN, como para el de LA.
searchGridSubCol <- expand.grid(gamma = c(0,0.01,0.1),
eta = c(0.1, 0.3, 0.5),
alpha = c(0.7,0.8,0.9))
set.seed(1234)
system.time(
rmseErrorsHyperparameters <- apply(searchGridSubCol, 1, function(parameterList){
#Extraemos los parámetros para testear
currentEta <- parameterList[["eta"]]
currentGamma <- parameterList[["gamma"]]
currentAlpha <- parameterList[["alpha"]]
param = list(eta = currentEta,
gamma = currentGamma,
alpha = currentAlpha,
eval_metric = "mae",
objective = "reg:linear",
booster = "gblinear")
#hacemos un modelo XGBoost Cross Validation
xgboostModelCV <- xgb.cv(data = dtrain_la, params = param, nrounds = 5000,
nfold = 5, showsd = TRUE,
metrics = "mae", verbose = T,
print_every_n = 1000, early_stopping_rounds = 100)
xvalidationScores <- as.data.frame(xgboostModelCV$evaluation_log)
mae <- tail(xvalidationScores$test_mae_mean, 1)
trmae <- tail(xvalidationScores$train_mae_mean,1)
stop_it <- tail(xvalidationScores$st)
output <- return(c(mae, trmae, currentEta, currentGamma, currentAlpha))}))## [1] train-mae:0.158863+0.000507 test-mae:0.158866+0.000834
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [430] train-mae:0.073647+0.000360 test-mae:0.073664+0.000995
##
## [1] train-mae:0.158862+0.000631 test-mae:0.158855+0.002740
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [415] train-mae:0.073643+0.000531 test-mae:0.073667+0.002089
##
## [1] train-mae:0.158862+0.000265 test-mae:0.158863+0.001473
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [428] train-mae:0.073646+0.000056 test-mae:0.073663+0.000311
##
## [1] train-mae:0.090079+0.000615 test-mae:0.090087+0.002502
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [141] train-mae:0.073627+0.000652 test-mae:0.073651+0.002451
##
## [1] train-mae:0.090078+0.000425 test-mae:0.090081+0.001815
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [140] train-mae:0.073624+0.000443 test-mae:0.073641+0.001711
##
## [1] train-mae:0.090078+0.000357 test-mae:0.090088+0.001181
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [142] train-mae:0.073627+0.000296 test-mae:0.073640+0.000854
##
## [1] train-mae:0.086735+0.000277 test-mae:0.086737+0.001074
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [84] train-mae:0.073627+0.000276 test-mae:0.073641+0.000914
##
## [1] train-mae:0.086735+0.000428 test-mae:0.086746+0.001759
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [85] train-mae:0.073625+0.000456 test-mae:0.073639+0.001705
##
## [1] train-mae:0.086735+0.000505 test-mae:0.086746+0.001677
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [86] train-mae:0.073628+0.000497 test-mae:0.073644+0.001516
##
## [1] train-mae:0.161137+0.000370 test-mae:0.161138+0.000743
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [440] train-mae:0.073651+0.000476 test-mae:0.073677+0.001503
##
## [1] train-mae:0.161137+0.000746 test-mae:0.161142+0.000957
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [427] train-mae:0.073648+0.000508 test-mae:0.073661+0.001645
##
## [1] train-mae:0.161137+0.000201 test-mae:0.161138+0.000832
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [430] train-mae:0.073644+0.000119 test-mae:0.073656+0.000692
##
## [1] train-mae:0.089828+0.000527 test-mae:0.089839+0.002192
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [138] train-mae:0.073625+0.000554 test-mae:0.073635+0.002182
##
## [1] train-mae:0.089827+0.000388 test-mae:0.089829+0.001440
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [144] train-mae:0.073625+0.000392 test-mae:0.073642+0.001461
##
## [1] train-mae:0.089827+0.000326 test-mae:0.089825+0.001372
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [140] train-mae:0.073630+0.000278 test-mae:0.073643+0.001128
##
## [1] train-mae:0.086736+0.000372 test-mae:0.086730+0.001527
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [83] train-mae:0.073624+0.000408 test-mae:0.073626+0.001486
##
## [1] train-mae:0.086736+0.000314 test-mae:0.086737+0.001307
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [83] train-mae:0.073629+0.000351 test-mae:0.073633+0.001317
##
## [1] train-mae:0.086735+0.000464 test-mae:0.086743+0.001848
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [84] train-mae:0.073625+0.000357 test-mae:0.073653+0.001394
##
## [1] train-mae:0.162782+0.000484 test-mae:0.162784+0.001763
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [427] train-mae:0.073645+0.000425 test-mae:0.073656+0.001385
##
## [1] train-mae:0.162783+0.000476 test-mae:0.162793+0.001464
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [437] train-mae:0.073647+0.000429 test-mae:0.073669+0.001525
##
## [1] train-mae:0.162783+0.000447 test-mae:0.162778+0.001771
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [438] train-mae:0.073645+0.000494 test-mae:0.073671+0.001665
##
## [1] train-mae:0.089682+0.000335 test-mae:0.089690+0.001085
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [139] train-mae:0.073630+0.000329 test-mae:0.073656+0.001023
##
## [1] train-mae:0.089682+0.000336 test-mae:0.089685+0.001188
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [141] train-mae:0.073627+0.000329 test-mae:0.073645+0.001134
##
## [1] train-mae:0.089682+0.000300 test-mae:0.089689+0.000915
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [142] train-mae:0.073630+0.000336 test-mae:0.073653+0.001373
##
## [1] train-mae:0.086781+0.000518 test-mae:0.086796+0.001956
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [82] train-mae:0.073622+0.000514 test-mae:0.073654+0.001857
##
## [1] train-mae:0.086779+0.000577 test-mae:0.086781+0.002195
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [83] train-mae:0.073625+0.000471 test-mae:0.073642+0.001673
##
## [1] train-mae:0.086779+0.000241 test-mae:0.086785+0.001068
## Multiple eval metrics are present. Will use test_mae for early stopping.
## Will train until test_mae hasn't improved in 100 rounds.
##
## Stopping. Best iteration:
## [83] train-mae:0.073624+0.000303 test-mae:0.073635+0.001197## user system elapsed
## 513.46 90.92 298.22Y vemos los resultados ordenados de forma creciente por TestMAE para ver los metaparámetros que dieron un MAE más bajo en el test.
set.seed(1234)
output <- as.data.frame(t(rmseErrorsHyperparameters))
varnames <- c("TestMAE", "TrainMAE", "Eta", "Gamma", "Alpha")
names(output) <- varnames
output[order(output$TestMAE),]## TestMAE TrainMAE Eta Gamma Alpha
## 19 0.0736604 0.0736490 0.1 0.00 0.9
## 12 0.0736616 0.0736490 0.1 0.10 0.8
## 13 0.0736632 0.0736524 0.3 0.00 0.8
## 11 0.0736672 0.0736530 0.1 0.01 0.8
## 3 0.0736704 0.0736520 0.1 0.10 0.7
## 6 0.0736706 0.0736582 0.3 0.10 0.7
## 1 0.0736708 0.0736532 0.1 0.00 0.7
## 16 0.0736712 0.0736686 0.5 0.00 0.8
## 15 0.0736716 0.0736578 0.3 0.10 0.8
## 5 0.0736720 0.0736544 0.3 0.01 0.7
## 2 0.0736722 0.0736474 0.1 0.01 0.7
## 14 0.0736722 0.0736546 0.3 0.01 0.8
## 20 0.0736732 0.0736520 0.1 0.01 0.9
## 23 0.0736742 0.0736530 0.3 0.01 0.9
## 21 0.0736758 0.0736498 0.1 0.10 0.9
## 27 0.0736774 0.0736662 0.5 0.10 0.9
## 17 0.0736798 0.0736738 0.5 0.01 0.8
## 24 0.0736798 0.0736572 0.3 0.10 0.9
## 4 0.0736822 0.0736576 0.3 0.00 0.7
## 22 0.0736830 0.0736546 0.3 0.00 0.9
## 10 0.0736838 0.0736574 0.1 0.00 0.8
## 26 0.0736838 0.0736656 0.5 0.01 0.9
## 8 0.0736890 0.0736746 0.5 0.01 0.7
## 7 0.0736914 0.0736756 0.5 0.00 0.7
## 9 0.0736958 0.0736786 0.5 0.10 0.7
## 25 0.0736974 0.0736648 0.5 0.00 0.9
## 18 0.0736976 0.0736694 0.5 0.10 0.8Vamos a tratar de hacer un modelo de XGBoost sin la validación cruzada, aplicándole los metaparámetros optimizados, y comparar los resultados con los modelos anteriores, para ver si la validación cruzada sobreajusta en exceso. Para Orange y Ventura primero.
set.seed(1234)
param <- list(booster = "gblinear",
objective = "reg:linear",
eval_metric = "rmse",
eta=0.1,
gamma = 0.00,
alpha = 0.7)
modelo_orve1 <- xgboost(params = param, data = dtrain_orve, nrounds = 5000, print_every_n = 1000, early_stop_round = 10, maximize = F, verbose = F)Ahora hacemos la predicción del dataset Test con nuestro modelo entrenado, para ver los resultados de las métricas.
prediccion1 <- predict (modelo_orve1,dtest_orve)
rmse_test_xgb <- RMSE( prediccion1, res_orve$test$log_error)
paste("Error de test (RMSE) del modelo XGB:", rmse_test_xgb)## [1] "Error de test (RMSE) del modelo XGB: 0.14772661688401"mse_test_xgb <- MSE(prediccion1, res_orve$test$log_error)
paste("Error de test (MSE) del modelo XGB:", mse_test_xgb)## [1] "Error de test (MSE) del modelo XGB: 0.021823153335995"MAE_test_xgb <- MAE(prediccion1, res_orve$test$log_error)
paste("Error de test (MAE) del modelo XGB:", MAE_test_xgb)## [1] "Error de test (MAE) del modelo XGB: 0.0558598835514259"MedianAE_test_xgb <- MedianAE(prediccion1, res_orve$test$log_error)
paste("Error de test (MedianAE) del modelo XGB:", MedianAE_test_xgb)## [1] "Error de test (MedianAE) del modelo XGB: 0.0285193836779952"Y ahora hacemos lo mismo para Los Angeles.
set.seed(1234)
param <- list(booster = "gblinear",
objective = "reg:linear",
eval_metric = "rmse",
eta=0.1,
gamma = 0.00,
alpha = 0.7)
modelo_la1 <- xgboost(params = param, data = dtrain_la, nrounds = 5000, print_every_n = 1000, early_stop_round = 10, maximize = F, verbose = F)Y de nuevo predecimos con el modelo la partición test de Los Angeles, para ver los resultados de las métricas.
prediccion2 <- predict (modelo_la1,dtest_la)
rmse_test_xgb <- RMSE( prediccion2, res_la$test$log_error)
paste("Error de test (RMSE) del modelo XGB:", rmse_test_xgb)## [1] "Error de test (RMSE) del modelo XGB: 0.15356373609137"mse_test_xgb <- MSE(prediccion2, res_la$test$log_error)
paste("Error de test (MSE) del modelo XGB:", mse_test_xgb)## [1] "Error de test (MSE) del modelo XGB: 0.02358182104234"MAE_test_xgb <- MAE(prediccion2, res_la$test$log_error)
paste("Error de test (MAE) del modelo XGB:", MAE_test_xgb)## [1] "Error de test (MAE) del modelo XGB: 0.0732008891672373"MedianAE_test_xgb <- MedianAE(prediccion2, res_la$test$log_error)
paste("Error de test (MedianAE) del modelo XGB:", MedianAE_test_xgb)## [1] "Error de test (MedianAE) del modelo XGB: 0.0370579702347517"importance_la <- xgb.importance(feature_names = train_sparse_la@Dimnames[[2]],
model = modelo_la1)
importance_orve <- xgb.importance(feature_names = train_sparse_orve@Dimnames[[2]],
model = modelo_orve1)
importance_orve## Feature Weight
## 1: calculated_finished_living_area 8.43821e-06
## 2: tax_amount -2.80833e-06
## 3: id_zip -1.46998e-06
## 4: lot_size_square_feet 1.28669e-07
## 5: id_city 9.96912e-08
## ---
## 292: week_transaction.49 0.00000e+00
## 293: week_transaction.50 0.00000e+00
## 294: week_transaction.51 0.00000e+00
## 295: week_transaction.52 0.00000e+00
## 296: week_transaction.53 0.00000e+00## Feature Weight
## 1: tax_amount -1.01109e-05
## 2: calculated_finished_living_area 9.89660e-06
## 3: id_county -5.29601e-06
## 4: tract_number -1.03016e-06
## 5: id_zip -9.48865e-07
## ---
## 292: week_transaction.49 0.00000e+00
## 293: week_transaction.50 0.00000e+00
## 294: week_transaction.51 0.00000e+00
## 295: week_transaction.52 0.00000e+00
## 296: week_transaction.53 0.00000e+00
5. Conclusiones
XGBoost obtiene un RMSE para los 2 condados de media de 0,14 y GBM de 0,16
XGBoost obtiene un MAE para los 2 condados de media de 0,06 y GBM de 0,065
Como era previsible XGBoost obtiene mejores resultados que GBM pero también hay que remarcar que los resultados no son muy difernetes entre sí, sin embargo el nivel de abstracción entre un modelo u otro si que puede ser relevante.